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Two-dimensional  models  of  three-component  Stirling  microrefrigerators  are 
examined.  The  models  consist  of  heater  and  cooler  ended  with  pistons  oscillating  out  of 
phase  and  connected  by  narrow  channels  with  a  linear  wall-temperature  distribution, 
simulating  the  regenerator.  Navier-Stokes  equations  governing  the  variable-density, 
oscillatory  flow  and  heat  transfer  within  two-dimensional  domains  with  moving 
boundaries  are  numerically  solved  by  using  the  SIMPLE  algorithm  incorporated  into  a 
moving  grid  technique  and  a  pressure-splitting  technique  that  is  newly  developed  in  this 
study.  Such  a  solution  methodology  is  accurate  and  efficient  for  the  complex  flow  and 
heat  transfer  problems  in  Stirling-type  devices.  The  pressure-splitting  technique  separates 
the  hydrodynamic  and  thermodynamic  pressures,  which  are  different  by  orders  of 
magnitude  for  gaseous  flows  with  heat  transfer  like  those  in  reciprocating  machines. 


Separation  of  the  pressures  increases  the  computing  accuracy  and  efficiency  substantially 
(with  saving  on  memory  by  50%  and  on  CPU  time  by  30%),  and  the  technique  can  be 
easily  incorporated  to  any  pressure-based  numerical  algorithms.  The  model's  cooling 
performance  is  sensitive  to  a  dimensionless  parameter  that  is  the  product  of  Prandtl 
number,  Womersley  number  squared,  and  two  dimensionless  geometric  parameters.  This 
parameter  should  be  kept  small  enough  for  the  regenerator  to  function  well.  As  the 
refrigeration  temperature  decreases,  the  cooling  capacity  decreases  linearly,  but  the 
power  input  increases  linearly,  and  consequently  the  performance  degrades.  It  is  found 
that  for  a  set  of  given  operation  conditions,  there  exists  a  lowest  temperature  below 
which  no  cooling  becomes  possible.  This  is  caused  by  the  increasing  cooling  loss 
through  the  channel  as  the  refrigeration  temperature  decreases.  Lengthening  the  channels 
decreases  such  loss,  but  meanwhile  increases  the  dead  space,  which  in  turn  reduces  the 
cooling  made  by  the  expansion.  The  net  effect  depends  on  these  two  competing  factors 
under  specific  conditions.  The  phase  angle  affects  the  cooling  capacity  and  efficiency, 
and  the  optimum  phase  angle  for  either  cooling  capacity  or  efficiency  is  near  90°.  For 
the  multichannel  model,  the  ratio  of  the  chamber  width  to  the  channel  width  also  affects 
the  cooling  performance,  and  the  optimum  ratio  depends  on  specific  operation  conditions. 


CHAPTER  1 
INTRODUCTION 


1.1    Background  of  Stirling  Refrigerators 
1.1.1    The  Stirling  Principle 

With  the  ban  on  the  production  of  chlorofluorocarbon  (CFC)  refrigerants  in  the 
U.S.  after  January  1,  1996  (and  eventually  the  ban  on  the  use  of  them  worldwide), 
alternative  cooling  devices  such  as  Stirling  refrigerators,  pulse  tube  refrigerators  and 
thermoacoustic  refrigerators  have  received  more  and  more  attention  (Chase  1995).  The 
Stirling  refrigerator  is  named  after  Scottish  clergyman  Robert  Stirling  who  invented  the 
Stirling  engine  in  1816.  A  Stirling-type  device  uses  a  common  gas  (not  ozone-depleting 
CFCs),  usually  air  or  helium,  as  the  working  medium  and  works  in  a  closed  cycle.  By 
compressing  and  expanding  the  gas  in  two  separated  chambers  at  different  temperatures, 
either  refrigeration  is  obtained  by  doing  work  on  the  gas  (as  in  a  refrigerator),  or 
mechanical  power  is  produced  by  consuming  heat  (as  in  an  engine).  Figure  1.1 
illustrates  the  working  principle  of  a  Stirling  cooling  device.  The  device  consists  of  three 
basic  heat-exchanging  components:  the  cold  chamber,  the  hot  chamber  and  the 
regenerator,  so  it  is  called  a  three-component  configuration.  As  shown  in  the  figure,  the 
chambers  are  connected  by  a  regenerator.  In  each  chamber,  there  is  a  piston  oscillating 
out  of  phase  with  the  piston  in  the  other  chamber.   The  regenerator  is  made  of  metallic 
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Fig.  1.1    Stirling  cycle  and  a  three-component  Stirling  refrigerator  configuration 


3 
porous  materials,  so  it  has  an  extremely  high  thermal  capacity  compared  with  the 
working  gas  and  meanwhile  has  excellent  thermal  contact  with  the  gas  in  it.  It  functions 
as  a  heat  sponge,  periodically  absorbing  and  releasing  heat  from  and  to  the  gas  flowing 
through  it.  The  Stirling  cycle  consists  of  four  idealized  processes  indicated  in  Fig.  1.1. 
Process  1-2-Isothermal  Compression: 

All  the  gas  is  located  in  the  hot  chamber  at  temperature  Th  and  the  right  piston 
compresses  the  gas.  As  a  result,  heat  is  rejected  from  the  gas  through  the  chamber  wall. 
Process  2-3-Constant  Volume  Precooling: 

The  two  pistons  move  in  unison  to  keep  the  total  volume  constant  and  to  move 
the  hot  gas  from  the  hot  chamber  to  the  cold  chamber  through  the  regenerator.    During 
this  process,  heat  is  stored  in  the  regenerator  and  the  gas  is  precooled. 
Process  3-4— Isothermal  Expansion: 

All  the  gas  is  in  the  cold  chamber  at  temperature  T„  and  the  left  piston  expands. 
Owing  to  the  expansion,  heat  is  absorbed  by  the  gas  through  the  chamber  wall.    Thus, 
refrigeration  is  achieved. 
Process  4-1— Constant  Volume  Preheating: 

The  pistons  move  in  unison  so  the  total  volume  is  constant  and  the  cold  gas  is 
moved  from  the  cold  chamber  to  the  hot  chamber  through  the  regenerator.  The 
previously  stored  heat  in  the  regenerator  is  now  released  back  to  the  cold  gas.  Then  the 
closed  system  is  ready  for  the  next  cycle. 

The  above  discussion  is  based  on  highly  idealized  conditions  under  which  the 
piston  motions  are  discontinuous  as  indicated  in  Fig.  1.1,  and  all  the  gas  is  in  the  hot 
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chamber  during  the  compression  and  all  is  in  the  cold  chamber  during  expansion.   The 

seidealized  conditions  also  imply  perfect  heat  transfer  and  infinite  heat  capacity  of  the 
regenerator  so  that  gas  temperature  in  either  end  chamber  is  maintained  constant  and  the 
compression  or  expansion  is  isothermal,  and  the  amount  of  heat  absorbed  by  the 
regenerator  is  exactly  the  same  as  that  restored.  It  is  not  so  in  reality.  First,  since  the 
discontinuous  piston  motions  are  mechanically  difficult  to  implement,  the  pistons  in  a  real 
Stirling  cooling  device  usually  oscillate  sinusoidally  such  that  the  hot-chamber  volume 
variation  leads  the  cold-chamber  volume  variation  by  a  phase  angle  of  about  90°. 
Second,  it  is  not  possible  for  all  the  gas  to  be  in  either  end  chamber  for  compression  or 
expansion  because  there  is  always  some  gas  in  the  regenerator  or  in  other  pockets  (the 
volume  that  can  not  be  swept  by  pistons  is  called  dead  space).  And  finally,  both  the 
thermal  capacity  of  the  regenerator  and  the  rate  of  heat  transfer  between  the  gas  and  the 
three  exchangers  are  finite.  Therefore,  the  gas  in  the  chambers  does  not  exactly  have 
the  wall  temperature  and  the  regenerator  may  not  function  perfectly  either.  As  a  result 
of  these  complicating  factors,  the  thermodynamic  cycle  is  not  precisely  composed  of  two 
isothermal  and  two  equal-volume  processes,  and  the  actual  P-V  diagram  will  be 
somewhat  different  from  the  one  shown  in  Fig.  1 . 1  (Chapter  4  will  show  a  realistic  P-V 
diagram). 

In  the  three-component  configuration,  the  expansion  or  compression  takes  place 
in  the  big  end  chambers,  so  there  may  not  be  enough  time  and  area  for  heat  exchange 
between  the  gas  and  the  chamber  walls.  For  better  performance,  usually  the  compression 
and  expansion  chambers  are  separated  from  the  heating  area  (heater)  and  the  cooling  area 
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Fig.  1.2  Schematic  of  a  five-component  Stirling  refrigerator 

(cooler),  respectively.  The  heater  or  cooler  is  made  up  of  many  metal  tubes  of  small 
diameters  such  that  the  heat  transfer  is  adequate.  Figure  1.2  demonstrates  a  five- 
component  configuration. 


1.1.2   Brief  History  of  Stirling  Machine  Development  and  Applications 

Details  on  the  development  history  of  Stirling-type  devices  can  be  found  in  books 
by  Walker  (1973),  Reader  and  Hooper  (1983),  West  (1986)  and  Organ  (1992).  Since 
Stirling's  invention  almost  two  centuries  ago,  many  different  Stirling  devices  have  been 
developed  and  found  their  various  potential  and  actual  applications:  automotive  engines, 
small  electric  generators,  underwater  power  units,  solar  thermal  conversion,  space 
power,  artificial  heart  power,  heat  pumping,  cryogenic  cooling  engines  (cryocoolers), 
near-ambient-temperature  refrigerating  machines  (refrigerators).  The  early  and  primary 
developments  focused  only  on  machines  providing  mechanical  driving  power  (i.e.,  prime 
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movers).  After  the  middle  of  the  nineteenth  century,  because  of  the  inventions  of  the 
internal  combustion  engine  and  electric  motor,  the  Stirling  prime  movers  gradually  lost 
their  commercial  market  to  these  new  inventions.  However,  the  appealing  characteristics 
of  the  Stirling  engine  still  attracted  engineers  to  continue  research  and  development. 
When  working  as  a  prime  mover,  a  Stirling  engine  can  use  essentially  all  kinds  of  heat 
sources:  solid,  liquid  or  gaseous  fossil  fuel  combustion,  solar  heat,  nuclear  heat,  etc. 
The  combustion  with  a  Stirling  engine  is  external  and  the  process  can  be  easily 
controlled,  producing  a  lower  air  pollution  level  than  an  internal  combustion  engine.  The 
Stirling  engine  does  not  require  any  valves  to  control  the  flow  or  pressure,  and  the 
combustion  (if  needed)  does  not  result  in  periodic  explosions.  So  the  engine  is  very  quiet 
and  the  machine  vibration  is  very  small  in  comparison  with  other  engines.  When 
working  as  a  refrigerator  a  Stirling  refrigerator  uses  no  CFCs,  so  it  is  environmentally 
safe. 

One  of  the  leading  institutions  for  the  Stirling  devices  is  the  Philips  Research 
Laboratories  in  Eindhoven,  Holland.  Back  to  the  late  1930s,  Philips  developed  portable 
Stirling  electric  generators  for  radio  sets.  Later  they  studied  and  applied  new 
technologies  to  the  then  existing  Stirling  engines.  By  the  mid-1940s,  they  had  increased 
the  Stirling  engine's  power  per  unit  weight  by  a  factor  of  50,  increased  its  size  per  unit 
power  by  a  factor  of  125  and  improved  its  efficiency  by  a  factor  of  15  (Reader  and 
Hooper  1983).    It  is  Philips  who  brought  the  Stirling  engine  into  its  modern  era. 

From  the  late  1950s,  Philips  along  with  several  European  and  American 
companies  (including  Ford  and  GM)  began  to  develop  Stirling  automotive  engines. 
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When  the  oil  crisis  hit  the  world  in  the  1970s,  the  interest  in  such  engines  increased 

suddenly  because  of  the  wide  variety  of  energy  (heat)  sources  that  could  be  used  by 
Stirling  engines.  In  1978,  sponsored  by  the  U.S.  Department  of  Energy,  several  U.S. 
companies  including  the  AM  General  Corporation,  Mechanical  Technology  Inc.  (MTI) 
and  the  United  Stirling  carried  out  intensive  researches  in  hopes  of  replacing  the  internal 
combustion  diesel  and  gasoline  engines  in  automobiles.  The  main  difficulty  with  Stirling 
automotive  engines  becoming  commercialized  is  that  they  cannot  compete  with  existing 
diesel  and  gasoline  engines  in  manufacture  and  maintenance  costs  because  the  latter  have 
already  been  well  developed  and  are  under  continuous  improvement.  As  a  result,  all 
Stirling  automotive  programs  stopped  about  a  decade  ago. 

One  of  the  active  areas  in  the  Stirling  research  and  development  is  the  artificial 
heart.  Since  the  1960s  the  U.S.  National  Institutes  of  Health  and  some  companies  (e.g., 
Stirling  Technology  Company)  have  been  developing  totally  implantable  artificial  heart 
devices  (West  1986,  White  1985).  Its  merits  of  quietness,  very  small  vibration  and  high 
reliability  make  the  Stirling  engine  the  best  power  provider  for  this  purpose.  A  small 
Stirling  engine  in  such  devices  uses  radioisotopes  as  the  heat  source  for  a  long  and 
continuous  life. 

Since  the  1960s  the  Stirling  engine  has  been  identified  as  the  key  to  the  future 
space  energy  conversion  systems  by  NASA  Lewis  Research  Center  (LeRC)  (Shaltens  and 
Schreiber  1989).  The  heat  sources  for  the  systems  can  be  nuclear  reactor  heat, 
radioisotopic  heat  or  solar  heat.  Stirling  systems  have  the  potential  of  meeting  space 
power's  requirements:  high  power  capacity,  high  thermal  efficiency,  high  reliability,  low 
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vibration  and  very  long  life.     NASA  initiated  the  free-piston  Stirling  space  power 

program  in  the  mid-1980s,  and  a  research  engine  was  built  by  MTI  under  contract  with 

NASA  (Tew  1988).     Some  main  design  goals  are:  output  electric  power  12.5  kW, 

thermal  efficiency  25%,  life  60,000  hours  (seven  years),  hot-end  temperature  1050  K, 

cold-end  temperature  525  K,  frequency  70  Hz  and  mean  pressure  150  atms.    Initial  test 

results  with  the  originally  built  engine  indicated  that  the  engine  performance  did  not 

reach  the  design  goal  (Tew  1988).  NASA  LeRC,  MTI  and  several  other  U.S.  companies 

have  done  intense  experimental  and  analytical  work  in  the  past  several  years  (Cairelli  et 

al.  1989-1991  &  1993,  Dochat  1990,  Jones  1990,  Dochat  and  Dhar  1991,  Thieme  and 

Swec  1992,  Wong  et  al.  1992),  and  some  significant  improvements  on  the  engine  have 

been  made  and  novel  component  technologies  have  been  obtained. 

The  commercialization  of  Stirling  solar  energy  conversion  systems  has  been 

pursued  for  almost  two  decades  in  several  countries,  e.g.,  Germany  (Keck  et  al.  1990), 

Japan  (Isshiki  etal.  1991),  Russian  (Loktinonov  1991),  Sweden  (Diver  et  al.  1990)  and 

the  U.S.  (Andraka  et  al.  1994).  In  such  systems  the  solar  heat  is  collected  by  a  parabolic 

dish  and  used  as  the  heat  source  of  a  Stirling  engine,  then  the  engine  can  drive  a  water 

pump  for  irrigation  or  drive  an  electrical  generator  for  electricity  supply.    The  Stirling 

solar  system  is  more  efficient  than  conversion  systems  using  photovoltaic  cells.    In  the 

early  solar  Stirling  engines,  collected  solar  rays  were  aimed  directly  at  the  engine  heater 

tubes.    This  produced  very  hot  spots  in  the  tubes  and  so  the  temperature  distribution 

there  was  not  uniform,  resulting  in  a  too  restrictive  requirement  on  materials  and  thus 

limiting  the  engine  life  and  performance.    To  solve  this  problem,  a  liquid-metal  reflux 
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receiver  is  used  to  transport  the  collected  solar  energy  to  the  heater  tubes.  In  the  reflux 
receiver,  a  liquid  metal  (usually  sodium  or  sodium-potassium  alloy  NaK-78)  evaporates 
at  the  solar  absorber  due  to  the  very  high  temperature  there,  and  then  condenses  at  the 
engine  heater  tubes.  The  latent  heat  released  by  the  condensing  vapour  is  utilized  to  heat 
the  tubes.  Most  of  the  current  technology  advances  are  made  in  the  U.S.,  primarily  due 
to  the  Department  of  Energy's  Solar  Thermal  Program  and  the  studies  made  by  some 
private  companies. 

The  most  successful  commercial  application  of  the  Stirling  engine  is  cooling. 
Seventeen  years  after  Stirling's  invention,  John  Herschel  in  1834  realized  that  the  Stirling 
machine  could  also  work  as  a  refrigerator,  and  the  first  Stirling  refrigerator  was  built  by 
Gorrie  in  1845  (Kohler  1960).  A  review  on  the  applications  of  Stirling  refrigerators  was 
given  by  Walker  et  al.  (1988).  In  1954  Philips  first  commercialized  a  Stirling  cryocooler 
with  a  cooling  capacity  of  1  kW  at  80  K  for  air  liquefaction.  Since  then  Stirling 
refrigerators  have  been  commercially  developed  for  various  applications  with  a  cooling 
capacity  from  a  few  milliwatts  to  several  kilowatts  and  with  refrigeration  temperatures 
from  as  low  as  a  few  degrees  Kelvin  to  near  ambient.  At  very  low  temperatures,  usually 
below  160  K,  the  efficiency  of  a  Stirling  refrigerator  is  much  higher  than  all  other 
cooling  devices  such  as  the  ones  using  the  Joule-Thomson  effect,  or  the  ones  using  the 
evaporation  of  a  liquid  (i.e.,  the  commonly  used  refrigerators).  So  Stirling  refrigerators 
have  dominated  the  world  market  of  gas  liquefiers  for  half  a  century.  Liquified  gases 
are  mainly  used  for  cryogenic  experiments,  in  which  the  temperature  is  usually  in  the 
cryogenic  range,  namely  below  120  K.     With  the  findings  of  more  and  more  high- 
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temperature  superconductors  and  their  increasing  applications  in  electronics,  the  current 

need  for  Stirling  cryocoolers  is  more  demanding  than  ever. 

Beside  the  large  cryocoolers  for  gas  liquefaction,  small  and  miniature  Stirling 
cryocoolers  have  also  been  widely  used  to  cool  some  special  small  sensors  and  electronic 
components  that  need  cryogenic  cooling  for  higher  resolution  and  sensitivity.  These 
include  infrared  detectors,  X-ray  and  7-ray  solid  state-detectors,  etc.  Increasing  numbers 
of  such  detectors  are  carried  by  airplanes,  space  shuttles,  space  stations,  satellites  and 
deep  space  probes  for  civilian  and  military  purposes  (e.g.,  Earth  observation, 
atmospheric  science,  astronomy  and  space  surveillance).  Since  the  1960s,  many  efforts 
have  been  made  in  the  U.S.  to  develop  spaceborne  cryogenic  systems  (Sherman  1982). 
The  required  cooling  power  for  the  systems  is  usually  between  a  few  milliwatts  and  tens 
of  watts  at  temperatures  below  80  K.  Due  to  the  restriction  on  the  size  and  weight  of 
the  loads  in  space  missions  (especially  the  long-term  ones),  it  is  not  feasible  to  carry  a 
large  amount  of  solid  or  liquid  cryogens  with  spacecraft.  Therefore  miniature 
mechanical  cryocoolers  (of  the  Stirling,  Joule-Thompson  absorption,  pulse  tube,  acoustic, 
turbo  Brayton  and  magnetic  types)  have  been  pursued  because  of  their  small  size  and 
lightweight,  mechanical  quietness,  and  expected  very  long  life  without  maintenance.  As 
reported  by  Chan  et  al.  (1990),  of  all  the  eight  known  space  flights  (between  1970  and 
1989)  using  mechanical  cryocoolers,  Stirling  cryocoolers  were  flown  six  times.  The  lead 
of  the  Stirling  cryocooler  over  others  is  attributed  to  the  mature  technologies  for  the 
large-size  Stirling  cryocoolers  as  mentioned  before  and  to  its  higher  efficiency  than  others 
in  the  temperature  range  of  interest. 
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In  search  of  alternatives  to  the  existing  CFC-containing  vapour-compression 
refrigerators  exclusively  used  at  near  ambient  temperatures,  the  Stirling  refrigerator  has 
recently  been  considered  as  a  promising  contender  (Walker  et  al.  1992,  Berchowitz  1992, 
Chase  1995).  According  to  Walker  et  al.,  at  near  ambient  temperatures,  the  coefficient 
of  performance  of  a  vapour-compression  refrigerator  is  about  60-70%  of  the  Carnot 
value,  whereas  that  of  the  Stirling  refrigerator  is  only  about  45%  of  the  Carnot  value, 
and  that  is  why  the  former  has  long  dominated  the  refrigeration  market  at  near  ambient 
temperatures.  Researchers  in  many  countries  are  now  developing  Stirling  refrigerators 
for  domestic  refrigeration,  air  conditioning,  electronic  cooling  and  so  on.  In  the  late 
1980s,  the  NASA  Manned  Spacecraft  Center  at  Houston,  Texas,  sponsored  a  company 
to  build  a  Stirling  refrigerator  for  manned  spacecraft.  Air  was  used  as  the  working  fluid 
in  the  device,  which  makes  it  safe  for  astronauts  in  case  of  refrigerant  leakage  and  simple 
for  gas  replenishment  since  the  working  gas  is  just  the  air  in  the  spacecraft.  To  replace 
the  existing  refrigerator  or  freezer  used  in  the  space  shuttle  mid-deck  for  the  preservation 
of  experiment  samples,  primarily  blood  and  urine,  NASA's  Life  Sciences  Projects 
Division  at  the  Johnson  Space  Center  in  1992  also  initiated  a  project  to  design  and 
develop  a  Stirling  refrigerator  (McDonald  et  al.  1994).  The  Stirling  refrigerator  was 
designed  by  Sunpower  Inc.  of  Athens,  Ohio.  This  unit  used  helium  as  the  working  fluid, 
and  was  lightweight  (101  lbs,  including  the  heat  pipes),  quiet  and  efficient  (heat  lift  of 
90.6  watts  at  4°C  with  a  power  input  of  60  watts,  and  heat  lift  of  88.5  watts  at  -22°C 
with  power  input  of  70  watts).  It  was  successfully  flown  on  board  a  space  shuttle  in 
February  1994  for  three  days  as  a  freezer  and  four  days  as  a  refrigerator.  Sunpower  Inc. 
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has  done  pioneer  work  in  developing  free-piston  Stirling  refrigerators  for  intermediate 
and  near  ambient  temperatures  and  in  commercializing  them  (Berchowitz  1992).  The 
company  has  developed  two  prototypes  of  Stirling  refrigerators:  "Super  Fridge"  and 
"Cool  Junior."  The  "Super  Fringe"  is  aimed  at  replacing  the  existing  domestic  vapour- 
compression  refrigerators,  and  is  capable  of  lifting  250  watts  at  -26°C  with  an  input  of 
215  watts.  The  "Cool  Junior"  is  very  small,  and  so  is  ideal  for  cooling  electronics, 
vaccines,  etc.  It  has  a  cooling  capacity  of  35  watts  at  -50°C  and  54  watts  at  -26°C. 
Work  on  the  development  of  Stirling  refrigerators  for  domestic  applications  has  also 
recently  been  done  by  Gold-Start  Co.  of  Korea  (Kim  et  al.  1993)  and  by  several  Japanese 
companies  and  institutions  like  Toshiba  (Otaka  et  al.  1993,  Isshiki  et  al.  1994).  Creare 
Inc.  of  Hanover,  New  Hampshire,  has  developed  a  proof-of-principle  diaphragm  defined 
Stirling  refrigerator  for  near  ambient  applications  (Shimko  et  al.  1994).  The  design 
innovation  in  this  machine  is  the  application  of  diaphragms  to  sweep  and  seal  the 
compression  and  expansion  volumes.  That  completely  avoids  the  serious  sealing  and 
lubricating  problems  associated  with  the  pistons  in  the  conventional  design. 

1.1.3   The  Challenges  in  Analyzing  and  Simulating  Stirling-Type  Devices 

In  reciprocating  machines,  such  as  internal  combustion  engines,  Stirling  type- 
devices  (prime  mover,  refrigerators  and  heat  pumps)  and  pulse  tube  refrigerators,  the 
flow  and  heat  transfer  are  time  periodic,  multidimensional  with  moving  boundaries  and 
with  abrupt  cross-section  changes.  Besides,  owing  to  large  temporal  volumetric 
oscillations  and  spatial  temperature  gradients,  density  varies  significantly  in  time  and  in 
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space,  and  pressure  level  also  oscillates  to  a  large  degree.   Moreover,  periodic  transition 

to  turbulence  may  also  occur  for  large  enough  machine  dimensions  or  at  high  enough 

operation  frequencies,  and  the  time-periodic  flow  and  heat  transfer  in  the  porous  medium 

of  the  regenerator  complicates  the  problem  even  further.     Therefore  analysing  or 

numerically  simulating  such  flow  and  heat  transfer  remains  a  challenging  task.    In  the 

design  and  analysis  of  such  machines,  usually  one-dimensional  (ID)  models  are  used 

along  with  empirical  formulas  for  gas-wall  heat  transfer  coefficients,  skin  friction 

coefficients,  pressure  drops  and  so  on  (Reader  and  Hooper  1983,  West  1986).    These 

formulas,  however,  are  correlated  only  from  the  experimental  data  of  incompressible 

(constant-density)  steady  unidirectional  flows;  and  it  has  been  shown  that  the  formulas 

and  consequently  the  ID  models  are  not  able  to  accurately  predict  the  performance  of 

Stirling  devices  (Geng  and  Tew  1992,  Smith  et  al.  1992,  Cairelli  et  al.  1993),  although 

they  meet  such  design  needs  as  simplicity  and  short  design  cycle. 

1.2  Literature  Review 
A  classical  thermodynamic  analysis  method  for  the  Stirling  engine  was  first 
presented  by  G.  Schmidt  in  1861  (Walker  1973).  The  differences  between  the  Schmidt 
theory  and  the  highly  idealized  Stirling  cycle  are  that  the  Schmidt  theory  allows  for  the 
more  realistic  sinusoidal  piston  motions,  and  that  it  also  takes  the  dead  space  into 
consideration.  This  theory,  however,  still  retains  some  of  the  highly  idealized  Stirling- 
cycle  assumptions.  For  example,  the  gas  temperature  distribution  in  each  heat  exchanger 
is  uniform,  therefore  the  thermodynamic  process  isothermal;  the  regeneration  is  perfect; 
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pressure  drops  across  engine  components  are  negligible  and  thus  the  instantaneous 
pressure  is  the  same  through  out  the  system.  Although  the  method  can  give  closed  form 
solutions  based  on  these  assumptions,  it  is  still  of  very  limited  use.  This  is  because  some 
assumptions  it  uses  are  still  far  from  the  real  situation,  and  so  the  predicated  engine 
characteristic  values  must  be  multiplied  by  a  factor  of  0.3  to  0.4  in  order  for  them  to 
match  those  of  a  real  engine. 

The  Schmidt  method  was  later  replaced  by  newer  analytical  methods  like  linear 
harmonic  analysis,  which  was  introduced  to  the  Stirling  community  after  1960s.  West 
(1984)  thoroughly  reviewed  the  linear  harmonic  analysis  method.  This  method  regards 
all  the  major  component-averaged  quantities  (volume,  pressure,  temperature,  etc.)  as 
unknown  sinusoidal  functions  that  are  linear  combinations  of  sin(wt)  and  cos(cot).  The 
governing  equations  used  in  the  method  are  algebraic  and  ordinary  differential  equations 
(ODEs)  written  for  each  individual  heat  exchanger  with  the  component-averaged 
quantities  as  the  dependent  variables.  Examples  of  the  algebraic  equations  are  the 
equation  of  state,  and  the  equation  of  continuity,  namely  the  summation  of  the  mass  in 
all  the  three  (sometimes  five)  components  equals  constant.  The  sinusoidal  functions  are 
substituted  into  the  governing  equations;  the  resulting  nonlinear  terms  of  sin(cjt)  and 
cos(wt)  are  linearized  by  using  binominal  expansion  first  and  then  simply  neglecting 
higher  order  terms.  Another  way  of  treating  the  nonlinearity  is  to  express  the  nonlinear 
terms  of  the  major  quantities  themselves  in  the  equations  as  Fourier  series  truncated  after 
the  sine  and  cosine  terms.  In  either  way,  the  linearized  equations  are  simply  a  few  (at 
most  less  than  dozens)  simultaneous  algebraic  equations.   The  latter  approach  is  able  to 
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account  for  the  effects  of  the  non- isothermal  or  transient  heat  transfer  process  and 

pressure  drops,  which  the  Schmidt  analysis  omits.  Huang  (1992)  reported  a  one- 
dimensional  computer  code  based  on  such  harmonic  analysis  method. 

In  the  past  decade,  fundamental  problems  of  oscillating  internal  flow  and  heat 
transfer  with  Stirling  applications  have  been  experimentally  studied  to  improve  analysis 
and  design  accuracy.  Tang  and  Cheng  (1993)  carried  out  experiments  on  incompressible 
periodically  reversing  pipe  flows,  and  correlated  a  cycle-averaged  Nusselt  number  in 
terms  of  the  Reynolds  number,  dynamic  Reynolds  number  and  the  nondimensional  fluid 
displacement.  Smith  and  co-workers  (Kormhauser  and  Smith  1989,  Smith  et  al.  1992) 
did  experiments  on  an  apparatus  including  a  piston-cylinder  space  connected  to  a 
dead-end  heat  exchanger  space,  and  proposed  a  "complex  Nusselt  number"  model  in 
which  heat  flux  consists  of  one  part  proportional  to  the  gas-wall  temperature  difference 
plus  a  second  part  proportional  to  the  rate  of  change  of  the  difference.  Tew  and  Geng 
(1992)  gave  a  review  on  related  experimental  projects  supported  by  NASA. 

Several  researchers  have  recently  begun  using  two-dimensional  (2D)  models  to 
simulate  Stirling  devices  numerically.  Gedeon  (1989)  developed  a  2D  computer  code 
modelling  the  unsteady  laminar  flow  and  heat  transfer  in  Stirling  components.  He  used 
Darcy's  law  to  simulate  the  flow  and  heat  transfer  in  the  regenerator,  and  applied  the 
density-based  Beam-Warming  scheme  (Beam  and  Warming  1978)  as  the  solution 
algorithm.  This  scheme  was  originally  designed  for  solving  compressible  (high  Mach 
number)  flows,  but  the  Mach  number  of  the  flows  in  a  Stirling  device  is  quite  low, 
typically  below  0.1.  (As  will  be  seen  later,  the  Mach  number  for  the  micro-configuration 
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in  the  present  work  is  below  1x10"*.)  Density-based  numerical  schemes,  e.g.,  the  Beam- 
Warming  schemes,  are  inefficient  and  inaccurate  in  solving  low  Mach  number  or 
incompressible  flows  (Karki  1986),  and  Gedeon  reported  similar  problems  encountered 
in  his  applications.  Ibrahim  and  co-workers  (1991,  1992)  also  used  2D  laminar  models 
to  study  Stirling  components  numerically,  and  suggested  a  way  to  correlate  the  gas-wall 
heat  transfer  coefficient  under  oscillating  flow  conditions,  and  they  investigated  the 
effects  of  sudden  cross-section  changes  on  heat  transfer  and  friction  coefficients.  The 
solution  method  employed  in  Ibrahim  and  co-workers'  work  was  the  pressure-based 
algorithm,  SIMPLE,  which  was  initially  developed  by  Patankar  (1980)  for  incompressible 
flows.  By  utilizing  a  2D  model,  Kurzweg  and  Zhao  (1993)  numerically  examined  the 
velocity  and  pressure  fields  only  within  one  end  chamber  of  a  Stirling  microrefrigerator 
with  various  orifice  configurations,  and  they  recommended  a  better  arrangement  of  the 
orifices.  They  also  used  the  SIMPLE  algorithm.  Density  in  the  chamber  was  regarded 
by  them  as  a  function  of  time  only. 

In  all  the  above-mentioned  2D  numerical  studies,  the  heat-exchanger  components 
were  separately  considered  as  open-ended  channels  (at  least  at  one  end),  and  the  velocity 
boundary  conditions  (B.C.)  at  the  component  entrance  or  exit  were  also  simply 
approximated  by  slug  flow  and  the  temperature  and  density  there  were  taken  as  uniform. 
Therefore,  the  effects  of  large  volume  and  pressure  oscillations  were  not  properly  taken 
into  account. 
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1.3  2D  Computational  Models  and  Scope  of  the  Present  Work 

The  main  purposes  of  this  work  are  1)  to  propose  2D  models  for  Stirling 
microrefrigerators,  2)  to  give  a  systematic  solution  methodology  suitable  for  solving  the 
variable-density  time-periodic  laminar  flow  and  heat  transfer  within  confined  regions  with 
moving  boundaries,  and  3)  to  numerically  investigate  the  flow  and  heat  transfer 
phenomena  in  such  configurations. 

For  computational  simplicity,  only  the  three-component  configuration  is 
considered  in  this  study.  The  model  adopted  here  integrates  all  the  three  heat  exchangers 
(cooler,  heater  and  regenerator)  together,  allowing  the  effect  of  large  pressure  and 
volume  variations  to  be  considered  properly.  To  reduce  the  computational  complexity, 
the  regenerator  is  replaced  by  narrow  channels  with  linear  wall  temperature  distributions. 
Since  the  channels  are  narrow,  there  is  enough  time  for  the  gas  flowing  in  them  to  adjust 
its  temperature  to  the  desirable  wall  temperature.  Two  types  of  models  are  studied: 
single-channel  models  Al  (Fig.  1.3  (a))  and  A2  (Fig.  1.3  (b));  and  multichannel  models 
Bl  (Fig.  1.4  (a))  and  B2  (Fig.  1.4  (b)).  In  the  single-channel  models  (Al  and  A2),  the 
two  chambers  are  connected  only  by  one  channel,  and  all  the  chamber  walls  are 
isothermal.  In  the  multichannel  models  (Bl  and  B2),  the  two  chambers  are  connected 
by  more  than  one  channel.  In  Models  A  and  B,  the  piston  surface  is  regarded  as 
insulating;  and  the  wall  temperature  of  the  left  chamber  is  set  at  Tc  and  that  of  the  right 
is  at  Th.  The  isothermal  surfaces  are  extended  into  part  of  the  connecting  channels  to 
a  length  "d"  in  Models  A2  and  B2.  As  indicated  in  the  figures,  the  channel  half  width 
is  "a",  the  total  channel  length  is  "L",  the  half  chamber  width  is  "b".    In  either  the  left 
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(a)  Single-channel  model:  A1 

-T.  T. 


X,e«  H 

xla(t  =  -L/2  -  o  +  esin(ut) 


x,ight  =  L/2  +  c  +  esin(cot+(po) 


(b)  Single-channel  model:  A2 


xlaft  =  -L/2  -  c  +  esin(wt) 


x,ight  =  L/2  +  c  +  esin(cot+(po) 


Fig.  1.3   Single-channel  models,    (a)  Model  Al;     (b)  Model  A2 
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or  the  right  chamber,  a  piston  oscillates  sinusoidally  with  an  amplitude  of  "fi"  and  an 
angular  frequency  of  "«".  The  piston  motion  in  the  right  (hot)  chamber  leads  that  in  the 
left  by  a  phase  angle  '>„"  for  cooling.  The  models  assume  a  small  clearance  of  "c-e" 
between  the  piston  and  the  vertical  wall.  The  reasons  for  this  are  1)  to  avoid  difficulties 
in  numerical  computations  if  there  is  no  clearance  and  2)  to  model  the  actual  situation 
where  there  is  always  a  clearance  and  some  dead  space.  Since  the  models  are  intended 
to  simulate  microrefrigerators,  the  models  are  in  the  centimeter  range. 

Even  though  the  flow  within  the  2D  models  is  laminar  because  of  the  very  small 
sizes  used,  numerically  solving  the  pertaining  heat  and  transfer  problem  still  constitutes 
a  difficult  task—the  challenging  elements  mentioned  before  remain  the  same.  Besides, 
in  the  variable-density  flows  with  heat  transfer,  there  are  significant  volumetric  changes. 
The  associated  thermodynamic  pressure,  therefore,  is  many  orders  larger  than  the 
hydrodynamic  pressure,  and  that  could  result  in  low  computational  accuracy  and 
efficiency,  and  even  physically  unrealistic  solutions.  To  overcome  this  problem,  a 
pressure-splitting  technique  is  developed  to  separate  the  thermodynamic  and 
hydrodynamic  pressures  in  this  study.  The  Navier-Stokes  and  energy  equations 
governing  the  variable-density,  oscillatory  flow  and  heat  transfer  are  solved  numerically 
by  using  the  pressure-based  algorithm,  SIMPLE  (Patankar  1980),  incorporated  with  the 
pressure-splitting  technique.  Since  the  domain  of  interest  is  bounded  by  moving  pistons, 
a  moving  grid  technique  as  used  by  Demirdzic  and  Peric  (1990)  is  also  incorporated  into 
the  algorithm.  This  study  proves  that  such  a  methodology  is  accurate  and  efficient  in 
solving  flow  and  heat  transfer  problems  like  those  occurring  in  Stirling  type  devices. 
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By  using  a  code  based  on  the  above  method,   the  2D  models  of  Stirling 

microrefrigerators  are  studied.  It  is  found  that  the  model's  regeneration  and  cooling 
performance  are  very  sensitive  to  a  dimensionless  parameter  that  is  the  product  of  the 
Womersley  number  squared,  Prandtl  number,  and  two  geometric  parameters. 
(Womersley  number  is  defined  as  a  =  a  yfulyfv,  where  01  is  the  oscillation  angular 
frequency  and  v  the  kinetic  viscosity.)  This  dimensionless  parameter  should  be  kept 
small  enough  to  let  the  gas  passing  through  the  channel  adequately  adjust  the  temperature 
to  that  of  the  wall's.  A  parametric  study  of  the  phase  angle  change  is  also  made  and  the 
optimum  phase  angle  is  shown  to  be  close  to  90°.  The  effects  of  varying  the  temperature 
ratios  between  expansion  and  compression  chambers  are  also  examined  in  this  study. 
The  cooling  capacity  is  found  to  vary  linearly  with  the  temperature  ratio.  Other  issues 
investigated  in  this  study  are  the  effects  of  changing  the  channel  length  and  changing  the 
ratio  of  chamber  width  to  channel  width.  Some  of  the  findings  from  the  current 
numerical  study  agree  well  with  experimental  data  with  actual  Stirling  refrigerators,  and 
that  suggests  the  proposed  2D  models  are  valuable  in  helping  gain  insights  into  the 
Stirling  microrefrigerators. 

The  arrangement  of  this  dissertation  is  as  follows.  The  governing  equations  and 
the  pressure-splitting  technique  will  be  discussed  in  Chapter  2.  Chapter  3  will  describe 
in  detail  the  numerical  techniques  solving  these  equations,  and  the  validation  of  the 
solution  methodology  will  also  be  made  in  Chapter  3.  Chapter  4  will  focus  on  the 
numerical  results  and  discussions,  and  finally  Chapter  5  contains  the  concluding  remarks. 


CHAPTER  2 
GOVERNING  EQUATIONS  AND  PRESSURE-SPLITTING  TECHNIQUE 


In  this  chapter,  the  governing  equations  of  fluid  mechanics  and  heat  transfer 
appropriate  to  the  problem  under  consideration  are  first  given.  A  method  is  devised  to 
split  the  pressure  into  hydrodynamic  and  thermodynamic  parts  so  that  the  computational 
efficiency  and  accuracy  for  incompressible  (very  low  Mach  number)  but  variable-density 
flows  with  heat  transfer  are  enhanced.  The  equations  and  their  boundary  conditions  are 
nondimensionalized.  Then  the  governing  equations  are  cast  into  a  form  that  is  suitable 
for  the  numerical  solution  algorithm  to  be  used.  At  the  end  of  this  chapter,  dimensional 
analysis  is  carried  out  for  the  cooling  power,  heating  power  and  the  rate  of  work  input 
in  the  Stirling  microrefrigerator  models. 

2.1  Governing  Equations 
The  governing  equations  of  viscous,  variable-density,  unsteady  flows  and  heat 
transfer  with  negligible  body  forces  in  domains  with  deformable  moving  boundaries  can 
be  written  in  an  integral  form  as  (Hansen  1967): 
Conservation  of  mass, 
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d 


V(t)  S(t) 

Conservation  of  momentum 
d 


—    f  pdV  +   i ptV-V^ndS  -  0  ,  (2.1) 

dt  A.  A, 


—   JpVdV  +  ^pVCV-V^-ndS  =  ^(-pI  +  f)-fldS  ,  (2.2) 

*  v(t)  s<t)  so) 


Conservation  of  energy, 


—  JpEdV+^pE(V-V^-fidS  =  -^a.-ndS+^V-(-pI  +  f)-ndS.      (2.3) 


d 
lv(t)  S(t)  S(t)  s<t) 


In  the  above  equations,  V(t)  represents  the  control  volume,  S(t)  the  moving  control 
surface,   y    tne  velocity  of  the  control  surface,  and  fi  the  outward  unit  normal  to  the 

control  surface;  p  is  the  density  of  the  fluid,   v  the  velocity  of  the  flow,  E  the  total 
energy  of  the  fluid  per  unit  mass,  p  the  pressure,  I   the  identity  stress  tensor,     i   the 

deviatoric  stress  tensor,  and  q  the  conduction  heat  flux  vector.    Here  the  assumption  of 
Newtonian  fluid  is  used. 

Notice  that  if  the  velocity  of  the  control  surface  is  the  same  as  the  fluid  velocity 
there,  then  the  Lagrangian  approach  is  taken;  if  the  velocity  of  the  control  surface  is 
zero,  then  the  Eulerian  approach  is  taken.  Otherwise  the  viewpoint  is  neither  Lagrangian 
nor  Eulerian,  as  happens  in  this  study. 

Since  the  body  force  is  negligible,  potential  energy  is  negligible,  too.  The  total 
energy  is  then  given  as 


E  =  e  +  -  V-V  ,  (2.4) 

2 
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where   i-V-V   is  the  kinetic  energy  per  unit  mass,  and  e  the  internal  energy  per  unit 
2 

mass.    The  internal  energy  for  a  perfect  gas  is  given  by, 

e  =  cvT.  (2.5) 

Here  cv  is  the  specific  heat  at  constant  volume,  and  T  the  absolute  temperature. 
According  to  Stokes'  hypothesis,  the  deviatoric  stress  tensor  is 

f  -  p.[VV+(VV)T]  -   |i-(V-V)I,  <26> 

where  the  superscript  "T"  represents  the  transposition  of  a  tensor,  and  n  is  the  molecular 
viscosity  of  the  fluid. 

The  equation  of  state  gives  the  relationship  among  p,  p  and  T,  and  reads 

p  =  pRT.  (2.7) 

where  R  is  the  gas  constant. 

In  the  energy  equation,  the  conduction  heat  flux  is  given  by  Fourier's  law 

q*  -  -k  VT  ,  (2-8) 

where  k  is  the  coefficient  of  thermal  conductivity. 

Viscosity  is  determined  by  using  Sutherland's  formula 

H  =  n    ZlIZe  (Z)3/2  (2.9) 

T  +  T       T 

in  which  Tr  is  some  reference  temperature,  fxr  the  viscosity  of  the  fluid  at  the  reference 

temperature,  and  T,,  a  constant  that  depends  on  the  fluid. 

Thermal  conductivity  k  can  be  determined  from  viscosity  /*  and  specific  heat  at 
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constant  pressure 

k  =  Ji^E  (2.10) 

Pr    ' 

where  Pr  is  the  Prandtl  number,  which  is  usually  a  constant  for  a  specific  fluid.    For  a 

perfect  gas, 

cp  -  cv  =  R,  (2.11) 

and  cp  =  7C„ ,  (2.12) 

where  y  is  called  the  ratio  of  specific  heats,  and  is  a  constant  for  a  given  gas. 


2.2  Separation  of  Hvdrodvnamic  and  Thermodynamic  Pressures 
The  pressure  p(x,y,z,t)  consists  of  hydrodynamic  pressure,  denoted  by  p\  and 
thermodynamic  pressure,  denoted  by  p',  i.e., 

p(x,y,z,t)  =  ph  +  p1.  (2.13) 

The  hydrodynamic  pressure  is  related  to  the  inertial  and  viscous  forces  and  is  usually  a 
function  of  both  time  and  space,  so 

ph  =  ph(x,y,z,t);  (2.14) 

and  the  thermodynamic  pressure  is  related  to  the  volumetric  changes  and  temperature 
variations.  At  very  low  Mach  numbers,  the  thermodynamic  pressure  is  a  function  of 
time  only,    i.e., 

P1  =  Pl(t),  (2.15) 

Consequently, 
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£^  -  0(M2)  «  1  .  (2.16) 

P' 

In  numerically  solving  the 

momentum      equation      (2.2), 

pressure  gradient  terms  like  dp/dx 

need  to  be  approximated.     As  a  P,  Pki 

1 1 1 —  x 


Fig.  2.1  Illustrative  ID  discretization 


demonstration,     consider    a     ID 

problem      where     the     pressure 

gradient  is  dp/dx.     As  shown  in 

Fig.  2.1,  the  pressure  gradient  at  xl+„,  =  (x<  +  xi+l  )/2,  denoted  as  (dp/dx)i+1/2,  is 

estimated  by 

(?W  -  ^^  •  (217) 

dx  x1+1-X[ 

where 

Pi  =  Phi  +  Pl,  (2-18) 

Pi+.  =  Phi+1  +  P',  (2.1?) 

From  the  above  discussion,  p'  is  independent  of  x  and  is  several  orders  larger  than  ph. 
So  what  is  done  in  Eq.  (2.17)  is  subtracting  two  very  close  numbers  from  one  another 
with  the  difference  between  them  being  many  orders  smaller  than  the  two  numbers 
themselves.  If  the  computer  used  does  not  have  enough  digits,  large  round-off  errors  and 
even  catastrophic  cancellation  may  occur  in  the  calculation  and  the  numerical  solutions 
are  erroneous  (Golub  and  Ortega  1992).    Although  using  double-precision  arithmetic  in 
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computations  can  increase  the  accuracy,  a  price  of  more  computer  memory,  storage  and 

CPU  time  must  be  paid. 

Wei  Shyy  *  suggested  that  effort  should  be  made  in  order  to  find  a  way  of 
separating  hydrodynamic  and  thermodynamic  pressures  so  that  the  computing  accuracy 
and  efficiency  could  be  increased.  The  rationale  follows.  Substituting  Eqs.  (2.18)  and 
(2.19)  into  Eq.  (2.17)  yields 

,dp,  (Pjti+p')-(Pih+p')       PiVPi"        dp\  (2  20) 

V'*"2"  x^.-x,  '    x^-x,    ~(dxW' 

Then  the  problem  of  catastrophic  cancellation  does  not  exist  anymore,  and  therefore  the 
computing  accuracy  is  enhanced  without  using  double-precision  arithmetic.  The 
computer  CPU  time,  memory  as  well  as  storage  are  saved.  Now  the  question  is  how  to 
determine  ph  and  p1. 

Such  a  method  is  developed  in  the  present  study.  First  consider  how  to  determine 
the  hydrodynamic  pressure  ph.  It  can  be  seen  from  Eq.  (2.20)  that  the  pressure  p 
appearing  in  the  momentum  equation  (2.2)  can  now  be  replaced  by  the  hydrodynamic 
pressure  ph  and  therefore  ph  is  solved  by  using  the  pressure-based  SIMPLE  algorithm, 
which  will  be  discussed  in  Chapter  3.  It  should  be  noted  that  the  total  pressure  appearing 
in  the  energy  equation  (2.3)  has  nothing  to  do  with  the  pressures  splitting  procedures 
because  the  pressure  there  is  replaced  with  pRT  via  the  equation  of  state  (2.7). 

Next  consider  the  thermodynamic  pressure  p1.   It  can  be  determined  by  using  the 


'Personal  communication,  1994. 
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equation  of  state  and  the  total  mass  of  the  whole  system,  and  the  procedure  follows. 
From  Eqs.  (2.7)  and  (2.13),  the  density  is  given  by 

JP_  =   Ph  +  P'  (2.21) 

RT  RT 

The  total  mass  of  the  entire  system  is  known  at  any  time  instant,  and 

Total  Mass  -  /  p(x!y>z,t)dV  -  /  £^Ml£W  dV  .  (2.22) 

I)  v(t)      RWA« 

Therefore, 


p'(t)  • [Total  Mass  -    f  |-iV]  . 

f—  dV  v(t)K1 


(2.23) 


V(.)RT 


In  this  way,  the  hydrodynamic  and  thermodynamic  pressures,  which  are  of 
different  physical  meaning  and  of  different  order  of  magnitudes,  are  determined 
separately,  and  consequently  the  computing  accuracy  and/or  efficiency  is  increased. 

2.3  Boundary  Conditions 
The  symmetry  of  the  geometry  and  the  boundary  conditions  under  consideration 
makes  the  solution  also  symmetric  about  the  x-axis  (see  Figs.  1.3  and  1.4).  Numerical 
calculation  takes  the  advantage  of  the  symmetry  by  considering  only  half  the  domain  of 
interest  as  the  computational  domain  with  symmetry  boundary  conditions  being  applied 
at  the  symmetry  line,  as  a  result  the  computing  effort  needed  is  greatly  reduced.  For 
example,  the  upper  half  domain  of  interest  in  the  present  work  is  taken  as  the 
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computational  domain.  In  what  follows,  the  boundary  conditions  for  Model  B2  (Fig.  1.4 

(b))  are  listed. 

Velocity  boundary  conditions  are  as  follows: 
x  =  xkn ,  0  i  y  i  b,  left  piston  wall,  no-slip: 

u  =  uiefi  =  -eusin(ait),    v  =  0.  (2.24) 

"icii  <  x  <  -L/2,  y  =  b,  imaginary  boundary,  symmetry: 

du/3y  =0,    v=0. 
x  =  -L/2,  a<  y  <b,  vertical  wall  of  the  left  chamber,  no-slip: 

u  =  0,   v  =  0. 
-L/2  <  x    <  L/2,  y  =  a,  upper  channel  wall,  no-slip: 

u  =  0,    v  =  0. 
x  =  L/2,  a  <  y  4  b,  vertical  wall  of  the  right  chamber,  no-slip: 

u  =  0,   v  =  0. 
L/2  <  x  <  xrlghl ,  y  =  b,  imaginary  boundary,  symmetry: 

du/dy  =0,    v=0. 
x  =  xright ,  0  <  y  i  b,  right  piston,  no-slip: 

u  =  "right  =  -ewsin(ut+^0),    v  =  0.  (2.25) 

xleft  <  x  <  x^t, ,  y  =  0,  symmetry: 

du/dy  =0.    v=0. 
The  temperature  boundary  conditions  are  as  follows: 
x  =  xM ,  0  <  y  <  b:  dT/dx  =  0. 

xlefl  <  x  <  -L/2,  y=b:  dT/dy  =  0. 
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x  =  -L/2,  a<y<b:  T  =  Tc . 

-L/2  <x  <  -L/2+d,  y=a:  T  =  Tu . 

-L/2+d  ui  L/2-d,  y=a:  T  =  Tc  +  (T„  -  TJ  (x  +  L/2  -  d)/(L  -  2d). 

L/2-d  <  x  i  L/2,  y=a:  T  =  Th . 

x  =  L/2,  a  <  y  *  b:  T  =  Th . 

L/2  <  x  <  xrlghl,  y=b:  dT/dy  =  0. 

x  =  xright ,  0  <  y  i  b:  dT/3x  =  0. 

*m,  <  "  <  X*u ,  y  =0:  3T/3y  =  0. 

Here  the  boundary  conditions  for  velocity  and  temperature  only  include  two  types, 
Dirichlet  and  Neumann,  and  all  the  Neumann  boundary  conditions  here  have  zero 
gradients.  Around  the  entire  boundary,  the  velocity  component  normal  to  the  boundary 
is  always  known. 

Pressure  does  not  requires  explicit  boundary  conditions  here.  As  will  be 
discussed  in  the  next  chapter,  the  numerical  algorithm  adopted  in  this  study  does  not  need 
explicit  pressure  boundary  conditions  when  the  velocity  normal  to  the  boundary  is  given. 

2.4  Nondimensionalization  of  Governing  Equations  and  Boundary  Conditions 
In  the  analysis  or  numerical  computation  of  a  physical  problem,  the  governing 
equations  and  boundary  conditions  are  usually  nondimensionalized.  The  reasons  for  this 
are  twofold:  first,  by  doing  so  the  nondimensional  parameters  governing  the  physical 
problem  are  obtained;  second,  computation  convenience  is  also  achieved  by  using 
nondimensional  quantities.     The  governing  equations  and  boundary  conditions  are 
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nondimensionalized  as  follows. 
First  reference  quantities  are  properly  chosen: 
Reference  length:  Lr  =  a  (channel  half  width); 

Reference  time:  tr  =  Uu  (co  =  angular  frequency  of  the  piston  oscillation); 
Reference  velocity:  Ur  =  LA  =  acJ-  Tne  reference  velocity  is  defined  in  terms  of  the 

reference  length  and  reference  time  and  cannot  be  arbitrarily  chosen; 
Reference  temperature:  Tr .    It  should  be  chosen  such  that  it  is  as  close  as  possible  to 

the  temporal  and  spatial  average  temperature; 
Reference  density:  p, .    Like  Tr ,  the  reference  density  is  also  chosen  so  that  it  is  as 

close  as  possible  to  the  temporal  and  spatial  average  density; 
Reference  viscosity:  /nr  =  n(Tr). 
Then  the  nondimensional  quantities,  which  are  denoted  by  bars,  can  be  defined  as 


u'  ~  -  ±        f.l 


*i 

xi     xi 

"  Lr  "    a  ' 

tr 

cot 

'  f* 

e 

e  ■  — . 

S--L- 

u,2 

au  pr  Tr 


S  +  -V-V, 


2  Ur 


Special  care  should  be  taken  in  nondimensionalizing  the  hydrodynamic  and 
thermodynamic  pressures  because  of  their  inherent  different  scales.  The  corresponding 
nondimensional  quantities  are  defined  as 
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prur2 


P«  -  -*- 


—    f  pdV  +   i p(VVb)-ndS  =  0  ,  (2.26) 

at  jL  £» 


P,RTr 

For  convenience  from  now  on  the  bars  over  the  nondimensional  quantities  will  be 
dropped  with  the  understanding  that  nondimensional  quantities  are  used,  unless  otherwise 
mentioned. 

Finally,  the  above  nondimensional  quantities  are  substituted  into  the  governing 
equations  and  boundary  conditions.    The  nondimensional  equations  are: 
Conservation  of  mass, 

V(t)  S(I) 

Conservation  of  momentum, 

—    fpVdV+  £ pVCV-V^-fidS  =      ^phndS  +  -^x-ndS  ,        (2.27) 
dt  v(t)  s(t)  s(t)  a  S(t) 

Conservation  of  energy, 

—  fpTdV  +   fpT(V-VJ-ndS  =  — *-  <f  V(p.T)ndS 

dt  *  "  Pr /v^  * 

v(t)  S(t)  rro  set) 

-(y-l)^pTV-iidS  ♦  ?(Y~1)M    <fv-x-ndS  (2.28) 

S(t)  °  S(l) 

♦  Y(Y-l)M2[-^/ipV-VdV  *   ^|pV-V(V-Vb)-ndS]  , 
dt  v(t)  2  s(t)  2 

Equation  of  state, 
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p<  ♦  YM2ph  =  pT  .  <2-29) 

Here  one  should  notice  the  following  two  points:  1)  the  pressure  in  the  momentum 
equation  (2.27)  is  just  the  hydrodynamic  pressure  p\  and  2)  the  form  of  the 
nondimensional  equation  of  state  (2.29). 

The  nondimensional  parameters  in  these  equations  are: 
Prandtl  number:  Pr  =  (/icp  Ik),  =  v,  Ik,  , 

Ratio  of  specific  heats:  7  =  cp  /cv , 

Mach  number:  M  =  ao)/(>RTr)"2  =  Ur/(7RTr)"2 

=  (reference  speed)  /  (speed  of  sound  at  Tr),  (2.30) 

Womersley  number:  a  =  a  (p^l/i,)1'2  =  a  (uilv,)m.  (2.31) 

This  parameter  characterizes  oscillatory  flows,  and  is  also  known  to  be  Stokes  number. 
Instead  of  the  above  parameter,  another  parameter  is  also  used  in  the  literature.  The 
latter  is  defined  as 

Va  =  a2w/er  =  a2  =  LrUr/><r.  (2.32) 

It  is  called  the  Valensi  number  or  nondimensional  frequency  (Ahn  and  Ibrahim  1991), 
or  oscillating  Reynolds  number  (Tew  1987). 

For  simplicity,  the  nondimensional  velocity  and  temperature  boundary  conditions 
will  not  be  listed  here.  The  nondimensional  parameters  resulting  from  the  boundary 
conditions  are  b/a,  c/a,  d/a,  e/a,  L/a,  Tc  IT, ,  Th  IT,  and  (p0 . 

Since  T  is  used  as  the  dependent  variable  to  be  solved  for  in  the  energy  equation 
(2.28),  quantities  p,  e  and  consequently  E  are  all  expressed  in  terms  of  T  there.  The 
second  term,  including  the  minus  sign  "-",  on  the  right-hand  side  (RHS)  of  the  equation 


34 
is  the  rate  at  which  work  is  done  to  compress  the  gas  and  therefore  the  rate  at  which  the 
internal  energy  increases  owing  to  compression;  the  third  term  on  the  RHS  is  the  rate  of 
dissipation  of  mechanical  energy  due  to  viscosity;  and  the  last  terms  in  the  square 
brackets  on  the  RHS  are  related  to  kinetic  energy.  In  this  study  the  dissipation  term  and 
the  terms  related  to  kinetic  energy  are  neglected  because  they  all  include  the  factor  M2 
and  M  here  is  always  smaller  than  10".  Therefore,  the  energy  equation  is  simplified, 
without  loss  of  computational  accuracy,  to 


-|pTdV+^pT(V-Vb)-ndS  =  ^^V(nT)-fldS-(Y-l)^pTV-ndS. 

*V(t)  S(l)  PTa   S(t)  S(t) 


(2.33) 


2.5  General  Form  of  the  Governing  Equations 
Any  of  the  integral  equations  presented  in  the  preceding  section  can  be  cast  into 
the  following  general  form 


—  Jp(|>dV  +   ^p<j)(V-Vb)-ndS  =   f(TV$yfidS  +   /fdV,         (2.34) 


d 

'v(t)  S(t)  SCt)  V(t) 

where  T  is  the  diffusion  coefficient,  f  the  source  per  unit  volume.  Each  specific  integral 
equation  corresponds  to  a  specific  set  of  <j>,  V  and  f,  which  are  tabulated  in  Table  2.1  for 
two  dimensions. 

The  first  term  on  the  left-hand  side  (LHS)  of  Eq.  (2.34)  is  the  rate  of  change  of 
the  extensive  property  inside  the  control  volume,  that  corresponds  to  the  intensive 
property  <t>;  the  second  term  on  the  LHS  is  the  rate  of  efflux  across  the  control  surface. 
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and  called  the  convection  term.   The  first  term  on  the  RHS  of  the  equation  is  called  the 

rate  of  diffusion,  and  the  second  term  on  the  RHS  the  source  term.  Anything  that  cannot 

be  put  into  the  three  terms  mentioned  above  is  put  into  the  source  term,  as  is  the  case 

in  the  momentum  and  energy  equations. 


Table  2.1  Terms  in  the  governing  equation  (2.34) 


<t> 

r 

t 

Continuity 

1 

0 

0 

u-Momentum 

u 

JL 

a2 

3ph     1   3  ,    3u.     2   9  .    cV      3,    dv\ 
-  -£—  + (n — ) ((i — )  +  — (u. — ) 

3x      3  3x3x      3dxdy      dy     dx. 
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2.6  Dimensional  Analysis 
The  physical  quantities  of  much  interest  with  the  computational  models  are: 
cooling  power,  heating  power  and  rate  of  work  input,  the  definitions  of  which  will  be 
given  in  the  next  chapter.  Since  they  all  have  the  dimension  of  power,  they  are  generally 
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denoted  as  V  for  the  dimensional  analysis  purpose.  For  the  flow  and  heat  transfer 
problem  under  consideration,  there  are  four  basic  dimensions,  [length],  [time],  [mass] 
and  [temperature].  Here  "[  ]"  means  "the  dimension  of".  The  four  basic  quantities, 
which  include  all  these  basic  dimensions  and  at  the  same  time  never  make  a 
nondimensional  quantity  themselves,  are  chosen  as  L, ,  t, ,  p,  and  Tr .  They  were  used 
as  the  references  in  the  previous  nondimensionalization  of  the  governing  equations  and 
boundary  conditions.  The  quantity  V  is  nondimensionalized  by  dividing  it  by  Lr5tr3prTt0, 
i.e.,  ,o,aV.    From  the  n  theorem  (Fox  and  McDonald  1978),  we  have 

VI  (p,aV)  =  ^{Pr,  7,  M,  a,  *>„,  Tc/Tr,  Th/Tr,  b/a,  c/a,  d/a,  e/a,  L/a}.       (2.35) 


CHAPTER  3 
NUMERICAL  TECHNIQUES  AND  CODE  VALIDATION 


The  integral  equations  discussed  in  the  last  chapter  need  to  be  solved  numerically 
for  the  specified  boundary  conditions.  The  gas  flow  under  consideration  is  still 
incompressible  even  though  the  density  of  the  gas  changes  considerably  with  space  and 
time.  This  is  because  the  Mach  number  here  is  very  small,  and  therefore  the  density 
variation  is  not  caused  by  pressure  gradients,  but  by  temperature  gradients  and 
volumetric  changes.  So  the  flow  should  be  called  incompressible  variable-density  flow, 
and  the  SIMPLE  algorithm  (Patankar  1980),  which  was  originally  developed  for 
incompressible  flows  and  heat  transfer,  can  be  adopted  as  the  proper  solution  algorithm. 
To  solve  the  problem  with  moving  boundaries,  the  moving  grid  technique  by  Demirdzic 
and  Peric  (1990)  is  incorporated  into  the  SIMPLE  algorithm.  The  pressure-splitting 
method  proposed  in  the  last  chapter  and  the  code  developed  accordingly  are  validated 
with  some  cases  whose  solutions  are  known  beforehand  so  that  the  code  can  be  used  for 
production  runs. 

3.1  Numerical  Grid  Generation 
In  order  for  the  integral  equations  to  be  solved  numerically,  the  entire  domain  is 
discretized  with  grids  into  a  system  of  subdomains,  or  cells.    At  every  time  step,  new 
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grids  are  generated  according  to  the  oscillating-piston  locations.  Owing  to  the 
rectangular  nature  of  the  boundaries  of  the  models  under  consideration,  Cartesian 
rectangular  grids  are  used.  Figures  3.1  (a)  and  3.1  (b)  show  two  typical  grid 
distributions  at  time  instants  ut  =  t/4  and  it,  respectively.  The  domain  is  decomposed 
into  three  blocks  with  each  in  the  left  chamber,  the  channel  and  the  right  chamber.  The 
total  number  of  grid  points  in  each  direction  within  a  particular  block  is  kept  constant  and 
the  grids  move  and  deform  only  in  the  x-direction  but  not  in  the  y-direction. 

Consider  the  grid  distribution  in  the  x-direction.  Suppose  there  are  N  grid  points 
to  be  distributed  along  a  line  segment  parallel  to  the  x-axis,  starting  at  x=x,  and  ending 
at  x  =  xN  ,  with  the  first  and  last  grid  spacings  being  specified  as  Ax„  and  Axf , 
respectively  (Fig.  3.2).  Grid  points  are  distributed  according  to  the  hyperbolic  tangent 
distribution  given  by  Thompson  et  al.  (1985)  as  follows. 
First  define 


y^A^f 


(3.1) 


N-l 

and 


H 


(3.2) 


Then  solve  one  of  the  following  nonlinear  equations  for  the  parameter  &  by  using  the 
Newton-Raphson  method  (Hornbeck  1982): 
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(a) 


(b) 


Fig.  3.1    Two  typical  grids,    (a)  wt=ir/4;  (b)  ut=i. 
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Fig.  3.2    ID  grids 
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sinhS        1 


5  G 


±         for  G  <  1  ,  (3-3) 


sin5       J_ 
~S~       G 


The  i-th  grid  is  distributed  at 


for  G  >  1  .  (3-4) 


(xN-x,)  ,  (3.5) 


1        '       H  +  (1-H)g, 

where 

tanh[8(— --)] 

gi.I{1 N-l     2    }.  (3.6) 

2  .    .8 

tanh— 

2 

For  each  chamber  block,  either  xN  or  x,  varies  with  time  due  to  piston 
oscillations.  The  two  end  spacings  are  given  in  terms  of  the  instantaneous  average  grid 
spacing  at  each  time  instant  so  that  a  good  moving  grid  distribution  is  achieved.  The 
hyperbolic  tangent  distribution  is  also  used  for  generating  the  grid  in  the  channel  block 
with  each  end  spacing  being  the  same  as  the  respective  adjacent  end  spacing  of  the 
chamber  block,  resulting  in  a  smooth  grid  distribution  near  the  chamber  and  channel 
interfaces.    Similar  procedures  are  applied  to  generate  grids  along  the  y-direction. 

3.2  Discretization  of  the  Governing  Equations 
As  mentioned  earlier,  the  governing  equations  need  to  be  discretized  on  a  finite 
number  of  cells  (i.e.,  finite  volumes),  so  that  they  can  be  solved  numerically.    The 
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staggered  grid  layout  is  adopted  here  to  avoid  the  pressure-velocity  decoupling  (see 

Patankar  1980).  Figure  3.3  shows  that  the  control  volumes  for  scalar  quantities  p,  p  and 
T  are  the  same  and  the  quantities  are  located  at  the  centers  of  the  control  volumes.  A 
typical  scalar  control  volume  is  highlighted  in  the  figure.  Velocity  components  (u  and 
v)  are  located  at  the  interfaces  of  the  scalar  control  volumes.  As  a  result,  the  u  control 
volumes  shift  in  the  x-direction  by  half  cell  distance  from  the  scalar  control  volumes,  and 
similarly  v  control  volumes  also  shift  in  the  y-direction  by  half  cell  distance.  Figure  3.4 
shows  a  typical  u  control  volume  and  a  typical  v  control  volume.  Each  of  the  quantities 
indicated  there  is  regarded  as  the  average  value  over  the  relevant  control  volume.  It  can 
be  seen  from  the  above  discussion  that  for  the  2D  staggered  grid  layout,  there  are  three 
sets  of  control  volumes:  one  set  for  all  scalar  quantities,  and  the  other  two  for  the  two 
velocity  components. 

The  general  equation  (2.34)  is  taken  here  as  an  example  to  demonstrate  the 
procedures  for  discretizing  the  governing  equations.  Consider  a  cell  centered  at  point 
P  with  neighboring  cells  centered  at  points  N,  W,  S  and  E  (as  shown  in  Fig.  3.5).  The 
midpoints  of  the  interfaces  between  cell  P  and  each  of  its  neighboring  cells  are  denoted 
by  lower  case  letters  n,  w,  s  and  e.  As  mentioned  in  the  last  section,  the  grid  points 
move  in  the  x-direction.  So  Fig.  3.5  also  shows  the  grid  velocities:  (ug)e~at  the  east 
interface,  and  (ue)w~at  the  west  interface.  The  grid  velocities  are  evaluated  by  using  first 
order  finite  difference  based  on  the  grid  coordinates  at  current  and  last  time  levels. 
From  Eq.  (2.34),  the  integral  equation  for  the  general  quantity  4>  over  cell  P  is 
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Fig.  3.3   Control  volumes  for  scalar  quantities 
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Fig.  3.4   Staggered  u  and  v  control  volumes 
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—  j  p<|)dV  +    £  p<|>(V-Vb)-ndS  -  ^rv<t>-  AdS  +    /  CdV  ,         (3.7) 

dtVJt)  Sp(t)  S(t)  Vp(t) 


•V^l)  Sp(t)  s<« 

where  the  subscript  "P"  denotes  the  quantity  of  cell  P. 
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Fig.  3.5  Schematic  of  a  general  cell  for  <t>  at  point  P 

Consider  the  rate-of-change  term  first.   In  reference  to  Fig.  3.5,  this  term  can  be 
discretized  as 


|/p(t,dV 

Vp(t) 


Pp4>pvp  ~  Pp<l>pV? 

At 


(3.8) 


where  At  is  the  time  step  size,  VP  the  cell  volume,  and  the  superscript  "o"  represents  a 
quantity  at  the  last  time  level.    Here  only  a  2D  case  is  considered,  so 
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VP  =  AP  =  AXpiyP ,  (3.9) 


(3.10) 


and 

consequently. 

V0 

PP<|)pAp 

At 

0  i  0  .    o 

Pp<PpAP 

The  convection 

term  is  discretized 

as 

4  P<t>(v-vb)-Ads  =  Pe<t>e(«-VeAyp-Pw<Mu-VwAyp 

sjw  (3.11) 

+  Pn<t>nVnAxpPs<t>sVsAxP  ■ 

Define  mass  fluxes  through  the  cell  faces 

Fc  =  [p(u-ue  )L  AyP ,     Fw  =  [p(u-ug  )]w  AyP ,  (3.12) 

F.  =  (pv)n  Axp  ,     Fs  =  (pv),  AxP,  0.13) 

then 


/  P(|)(V-Vb)-fidS  -  Fe<),e  -  Fw<t.w  ♦  F,4s,  -  ?.♦.  .  (3.14) 

Spffl 

The  diffusion  term  is  then  discretized  as 


/rvcMdS  -  (r^)eAy]j-(r£)wAyp 

*#>  (3.15) 

+  (r^)Ax  -(r^t)Ax  . 
ay  "     p        3y  "     ' 

The  first  derivatives  in  the  above  expressions  are  approximated  by  first  order  finite 
differencing,  e.g., 
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a*    .  Ob-o,  (3 16) 

ax h        8x„ 


Thus  one  has, 


RHS  of  Eq.(3.15)  -  r.-^^-^-r.-^^-^) 

Ax„  Ax„  o  17) 

=  De(<|)E-<l>p)-Dw(<|)p-<t)w)+Dn(<t)N-4.p)-D3(<|)p-(|)s)  , 


where  the  diffusion  fluxes  are  defined  as 


Avp  Ayp  Axp  Axp  ,->  io> 

D.'r^'D*"r^-    D--r«^'    D'"r'^-       <318) 

The  discretization  of  the  source  term  is  straightforward, 


/  fdV  -  VPCP  -  APCP  .  (3.19) 

Vp(t) 

Substituting  all  the  discretized  terms  into  the  integral  equation  (3.7)  gives 
Ppap4.p-PpAp°4,p  +Fd>    F<t>   = 

At  e^e    *WW    'nfn     *  s^s  (3.20) 

De(<t»E -fy) -Dw(<|)p -<|>w) +Dn(<|>N  (|)p)-Ds(<t)p -<t>s)  +ApCp  . 


Notice  that  the  fully  implicit  first  order  backward  Euler  scheme  is  used.  In  the  above 
equation,  the  4>  values  at  cell  faces,  i.e.  4>c ,  4>„ ,  <£w  and  0S,  have  to  be  expressed  in 
terms  of  the  values  at  cell  centers,  i.e.  <t>E  ,  <£N  ,  <£w  and  4>s,  since  the  latter  are 
regarded  as  unknowns  to  be  solved.    Convection  schemes  are  needed  to  do  this  (for 
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details,  see  Patankar  1980,  Shyy  1994).    Applying  a  convection  scheme  to  the  above 

equation  and  then  collecting  terms  gives  the  algebraic  equation  for  4>  in  the  final  form: 
M>p  "  aE<t>E+aw(t,w+aN<t>N+as<t)s+B  "  £anb<t>nb  +  B  -  (3.21) 

nb 

where,  for  compactness,  the  subscript  "nb"  denotes  evaluations  at  the  neighboring  points 

of  point  P,  and  the  coefficients  are  defined  by 

aE  =  DeC(|Pece|)  +  max{-Fe,  0},  (3-22) 

aw  =  DwC(|Pecw|)  +  max{Fw,  0},  (3.23) 

aN  =  D„C(  |  Pec„ | )  +  max{-F„ ,  0} ,  (3.24) 

as  =  DsC(|Pecs|)  +  max{Fs,  0},  (325) 

B  =  CpAp  ♦  aP°(t)?  ,  0.26) 


o       (pA£ 

ap  - 


(3.27) 


At 
ap  =  aE+aw  +  aN+as+aP°  =  Eaab  +  aP  ■  (3-28) 

nb 

The  symbol  "Pec"  in  these  expressions  is  the  cell  Peclet  number,  which  is  the  ratio  of 
the  convection  flux  to  the  diffusion  flux.  The  Peclet  number  at  the  east  face,  Pece ,  is 
defined  as 

Pec6  =  Fe  /  De .  (3.29) 

The  function  C  depends  on  the  convection  scheme  used,  and  in  this  work  the  hybrid 
scheme  is  used.   Then  C  takes  the  form 

C(  |  Pec  | )  =  max{0,l-(V4)|Pec|}.  (3.30) 
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This  scheme  is  essentially  central  difference  for  |Pece|  <  2,  and  switches  to  first-order 

upwind  for  |Pece|   >  2.    For  the  problem  under  consideration,  the  source  term  is  only 

linear  in  4>  at  most,  and  so  for  simplicity  it  is  just  lumped  into  the  term  B  (Eq.  (3.26)). 

Special  care  should  be  taken  in  the  discretization  of  boundary  cells.  As  discussed 
before,  there  are  only  two  kinds  of  boundary  conditions  involved:  Dirichlet  and 
Neumann.  Figure  3.6  shows  cell  P  with  the  south  face  being  the  boundary  where  the 
Dirichlet  boundary  condition  is  specified  there,  namely,  <t>,  is  given.  In  this  case,  let 
point  S  coincide  with  point  s,  then 

<t>s=<t>,.  (3.31) 

Here  also  note  the  special  definition  of  <5y,  at  the  boundary  cell,  as  could  be  seen  clearly 
by  comparing  it  with  that  in  Fig.  3.5.  With  these  two  special  definitions,  the  procedures 
and  formulations  to  discretize  the  boundary  cell  are  exactly  the  same  as  those  for  an 
interior  cell  as  discussed  before.  By  the  same  token,  all  other  boundary  cells,  whether 
they  are  at  the  north  or  the  east  boundary,  are  discretized. 

The  Neumann  boundary  condition  is  treated  in  another  way.  Again,  a  cell  P  with 
the  south  face  being  the  boundary  where  the  d4>/dy  is  specified  is  taken  as  an  example 
to  demonstrate  the  boundary  treatment.  As  shown  in  Fig.  3.7,  an  imaginary  cell  S  is 
added  to  the  south  of  cell  P  such  that 

Ays=AyP,  (3.32) 

and  so  5ys=AyP.  (3.33) 

Applying  the  boundary  condition  with  central  difference  approximation  yields 
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Fig.  3.6  A  south  boundary  cell  with  Dirichlet  boundary  condition 
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Fig.  3.7   A  south  boundary  cell  with  Neumann  boundary  condition 
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34,      _  <Pp^s  (3.34) 

Vs  8ys     ' 


then. 


<Ds  =  4»p-(^)s8ys.  (3-35) 

dy 
Now  the  boundary  cell  P  can  be  treated  as  if  it  were  an  interior  cell.    As  mentioned  in 
the  last  chapter,  all  the  Neumann  boundary  conditions  are  with  zero  gradients  in  this 
study,  consequently  the  last  equation  gives  $,  =  <j>,  for  such  special  cases. 

Throughout  the  discretization  process,  linear  interpolation  is  frequently  used  to 
evaluate  values  whenever  they  are  not  readily  available,  except  <t>'s  themselves  at  cell 
faces,  which  are  determined  by  convection  schemes  instead.  For  instance,  in  evaluating 
the  mass  flux  F„  =  /onv„AxP  (refer  to  Fig.  3.7),  if  pa  is  not  readily  available,  as  is  the 
case  when  the  control  volume  is  for  any  quantity  but  v,  then  the  density  there  should  be 
approximated  by  linearly  interpolating  the  values  at  the  neighboring  points. 

It  has  been  seen  that  the  involvement  of  moving  grids  doesn't  complicate  the 
algorithm  much.  It  is  the  velocity  relative  to  the  control  surface  that  is  used  in  defining 
the  convection  mass  flux  (refer  to  Eqs.  (3. 12)  and  (3. 13)).  Since  the  grid  velocity  in  this 
work  is  neither  zero  nor  the  fluid  velocity,  the  grid  is  neither  Eulerian  nor  Lagrangian. 

3.3  Solution  Procedures 
With  the  SIMPLE  algorithm,  the  discretized  governing  equations  are  solved  one 
at  a  time  in  a  sequential  manner.     Besides,  pressure  is  used  as  the  basis,  i.e.,  the 
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hydrodynamic  pressure  is  determined  by  letting  the  velocity  field  satisfy  the  continuity 

equation,  as  is  in  contrast  to  the  compressible  algorithms  (also  called  density-based 

algorithms)  that  use  density  as  the  basis  by  solving  the  density  from  the  continuity 

equation  and  then  determine  pressure  from  the  equation  of  state. 

At  the  very  beginning  of  the  calculation,  initial  conditions  must  be  properly 
specified.  Then  time  is  increased  by  At;  at  the  new  time  level  grids  are  generated  based 
on  the  new  boundaries  and  boundary  conditions  are  specified  accordingly.  The  flow  and 
thermal  fields  need  to  be  calculated,  and  the  procedures  are  described  as  follows. 

First,  all  the  unknown  quantities  are  guessed,  and  denoted  as  u\  v*,  (p1)*,  (pb)*, 
T*.  p*,  etc.  Good  initial  guesses  should  correspond  to  the  values  at  the  last  time  level. 
With  these  guessed  quantities,  the  u-momentum  equation  is  discretized  first.  Since  the 
way  to  discretize  a  general  equation  was  already  given  in  detail  in  the  last  section,  only 
the  final  form  of  the  discretized  u-momentum  equation  for  an  u-control  volume  P  is  given 
here. 

<¥»p  =  £  «W»b  -  aXp(Pe  "Pp  )  +s„  .  (3.36) 

no 

where  the  source  term  is  intentionally  written  into  two  parts  such  that  the  part  solely 
related  to  the  driving  force~pressure~is  separated  from  all  other  contributions,  as  this 
will  help  the  discussion  of  pressure- velocity  coupling. 

Equation  (3.36)  is  linearized  already  because  the  coefficients  involved  are 
evaluated  in  terms  of  the  guessed  u*  values  but  not  the  current  u  values.  Every  cell 
results  in  an  equation  like  Eq.  (3.36),  and  thus  a  system  of  linear  algebraic  equations  for 
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u  is  formed,  which  is  solved  by  an  iterative  method-the  line  successive  relaxation 
method  (Hirsch  1988).    The  u*  values  can  now  be  updated  by  the  newly  obtained  u 
values, 

uf  =  w„u*  +  (l-wu)u,  (3.37) 

where  u„  is  an  underrelaxation  parameter  (0<a)„<  1).   Owing  to  the  nonlinear  coupling 
between  the  u-momentum  equation  and  other  equations,  underrelaxation  would  be  needed 
in  updating  u*  values  so  that  the  solution  procedure  converges. 
The  discretized  v-momentum  equation  is  given  as 

3pvp  =  £  Wat, "  Axp(Pn    Pp")  +sv  •  (3.38) 

nb 

In  a  similar  way,  v  is  solved  and  v*  is  updated  with  underrelaxation. 

The  velocity  field  determined  this  way,  however,  may  not  necessarily  satisfy  the 
continuity  equation  because  the  solution  is  based  on  a  guessed  hydrodynamic  pressure 
field  (ph)*.  Therefore,  the  discretizing  the  continuity  equation  over  a  pressure-control 
volume  at  point  P  (refer  to  Fig.  3.5)  yields 

(P*A)p-(pA)°  (339) 

+Fe  -Fw+Fn-Fs    -  S*  (»  0)  .  V-W 

At 

The  pressure  needs  to  be  corrected  so  that  the  resultant  velocity  field  satisfies  the 

continuity  equation.      If  the  pressure  correction   is  denoted  by  p',   then,   from  the 

momentum  equations  (3.36)  and  (3.38),  the  corresponding  velocity  corrections  at  the 
interfaces  of  the  pressure-control  volume  are  given  by 
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ue  ■  "^(Pb-Pp)  ,      u,  -  -^(Pp-Pw)  ,  (3-40) 

ap  ap 

and 


Ax„      .       .  .  Axp      . 

vn  -  "  — (Pn-Pp)  ,      vs  =  -_^(pp-ps)  .  (3.41) 

ap  ap 

Since  there  is  a  half-cell  shift  between  the  velocity-  and  pressure-control  volumes,  an 

index  shift  is  also  applied  when  deriving  the  above  equations  from  the  momentum 

equations.    The  coefficient  a,,  in  Eqs.  (3.40)  differ  from  those  in  (3.41)  because  the  u- 

and  v-control  volumes  are  not  collocated.    The  corresponding  flux  corrections  are 


Ayp      ,       , 

Fe    '    PeU'eAyp    "    "P. (PE"Pp)AyP  . 


(3.42) 


Ayp      .       . 

Fw  *  Pwu'wAyp =  -Pw — (Pp-Pw)AyP . 
ap 

Axp     , 

F«    =    PnV'nAxP    =    "Pn (Pn"Pp)Axp   . 

3P 

Axp      , 
FS  *  PsV'sAxP  "   "P. (Pp"Ps)Axp  • 

ap 


It  is  required  that  the  corrected  pressure  and  velocity  satisfy 

(p-A)p-(PA)g  ,      ,       .      ,  (343) 

.  *e       "w     *n      's       *e     *w     *n     *s         u    ' 

Substituting  Eqs.  (3.42)  into  Eq.  (3.43)  and  collecting  terms  give  the  pressure  correction 
equation 
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Vp    "    £bnbPnb-Sm  .  (3-44) 

Db 

where  b's  are  the  coefficients,  and  s*m  is  defined  in  Eq.  (3.39).    The  latter  quantity 
serves  as  a  measure  of  how  closely  the  continuity  equation  is  satisfied. 

It  was  mentioned  in  the  last  chapter  that  there  are  no  explicit  pressure  boundary 
conditions  for  the  problems  under  consideration.  Pressure  correction  equations  like  Eq. 
(3.44),  however,  do  need  boundary  conditions.  Boundary  conditions  for  p'  need  to  be 
derived  from  the  given  velocity  boundary  conditions.  For  example,  consider  a  cell  with 
the  west  face  as  the  boundary  where  the  normal  velocity  component,  u,  is  predefined, 
then  there  should  be  no  velocity  (and  therefore  mass  flux)  correction  there,  i.e.,  F„'  in 
Eqs.  (3.42)  is  zero.  The  Fw'  term  will  not  appear  in  Eq.  (3.43),  so  the  coefficient  bw 
in  Eq.  (3.44)  takes  the  value  of  zero  (bw  =  0). 

After  pressure  corrections  are  solved  from  the  pressure  correction  equation,  the 
hydrodynamic  pressure  is  updated  by 

ph  =  (ph)'  +  u„p\  (3.45) 

where  underrelaxation  is  used  again  to  maintain  convergence  (0<up<  1). 

The  corresponding  velocity  corrections  u'  and  v'  are  calculated  from  Eqs.  (3.40) 
and  (3.41),  respectively,  and  the  velocity  components  are  then  updated  as 

u  =  u"  +  u\  (3.46) 

v=v*+v'.  (3.47) 

One  should  note  that  the  energy  equation  is  coupled  with  the  continuity  and 
momentum  equations,  and  the  reasons  are  1)  viscosity  is  a  function  of  temperature,  and 
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2)  for  the  variable-density  problem,  density  is  also  a  function  of  temperature.  So  the 
discretized  energy  equation  is  solved  with  the  guessed  density  p"  and  the  above  updated 
velocity  components.  The  thermodynamic  pressure  p1  can  then  be  computed  by  using  the 
method  described  in  the  last  chapter  via  Eq.  (2.23),  and  density  is  also  updated  from  Eq. 
(2.21).    Underrelaxation  is  also  used  in  updating  density. 

With  all  the  updated  quantities,  one  goes  back  to  the  momentum  equation  again 
and  repeats  the  solution  procedures  until  the  solutions  are  fully  converged  at  the  given 
time  level.  The  time  level  is  next  increased  by  one  when  solutions  have  converged. 
Since  the  time-periodic  problem  is  being  simulated,  it  would  take  several  cycles  (or 
periods)  for  the  numerical  solution  to  reach  a  time-periodic  state.  Figure  3.8  gives  a 
summary  of  the  overall  solution  procedures  in  a  flowchart. 

3.4  Initial  and  Slow-Start  Boundary  Conditions 
Because  of  the  iterative  feature  of  the  solution  procedure,  good  initial  guesses  for 
the  unknown  quantities  can  make  the  solution  converge  fast,  as  is  especially  the  case  at 
the  very  first  time  step.  At  t=0,  the  left  piston  velocity  is  zero  (Eq.  (2.24)),  but  the 
right  piston  velocity  is  not  if  <p0  ^  0  (Eq.  (2.25)).  Meanwhile,  if  Tc^Th,  it  would  be 
difficult  to  specify  the  precise  initial  velocity  and  temperature  conditions  for  the  solution 
process  to  start  with  because  the  initial  velocity  and  temperature  fields  are  not  uniform 
in  such  conditions.  To  circumvent  such  a  difficulty,  a  slow-start  procedure  was  utilized 
in  specifying  temperature  boundary  conditions,  and  the  phase  angle  lag  <p,  and  therefore 
the  velocity  boundary  conditions  as  well. 
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Fig.  3.8  Solution  procedures  for  variable-density  time-periodic 
flow  and  heat  transfer  problems 
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A  quantity  <p  is  defined  by  the  following  expression, 


l2[l-cos(2ot)|  ,    0S«t<4 

<P      ,       —  <  Olt    , 


(3.48) 


where  p„  is  the  desirable  phase  angle.  The  newly  defined  ip  will  take  the  place  of  <pa  in 
the  expression  for  xr^h,  in  Figs.  1.3  and  1.4,  and  also  in  Eq.  (2.25)  for  u^,.  Then  at 
t=0,  both  pistons  have  zero  velocity,  and  so  the  velocity  in  the  whole  computational 
domain  is  also  zero.  It  takes  a  quarter  of  the  first  cycle  for  if  to  adjust  smoothly  to  its 
final  value  <p0 . 

For  temperature  boundary  conditions,  suppose  that  the  desirable  chamber  wall 
temperatures  are:  Tc  for  the  left  chamber  and  Th  for  the  right,  then  define 


^!(Th+T)-(Th^Tc)[  1-008(2^)]) 


0  <  ait  < 


<ait 


(3.49) 


and 


i{(Th+T)  +  Hilll^[l-cos(2wt)]}  ,    0<ort<l,      (35Q) 

T.  ,    -  <  art  , 

h  2 


where  Tlefl  and  Tm,  are  the  temperatures  to  be  specified  as  the  boundary  conditions  for 


the  corresponding  chamber  walls.      At  t=0,   Tleft=T^ 


hence  a   uniform   initial 
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temperature  field   is  expected.      Over  a  quarter  of  the  first  cycle,   the  boundary 

temperatures  gradually  develop  to  the  desirable  values. 

3.5  Code  Validation 
The  proposed  pressure-splitting  method  and  the  code  developed  were  validated  for 
variable-density,  time-periodic  flows  and  heat  transfer  in  regions  with  moving 
boundaries.  The  first  test  case  is  a  2D  piston-cylinder  set-up  with  all  the  boundaries 
being  insulated  (Fig.  3.9).  Air  is  present  in  the  cylinder  with  initial  pressure  p„=l  std 
atm  (  =  1.013xl03  N/m2),  temperature  T0=288.15  K,  and  density  p„=  1.225  kg/m3.  At 
the  very  beginning,  the  piston  stands  still  at  x=0,  so  there  is  no  fluid  motion.  Then  the 
piston  oscillates,  and  its  location  is  given  by 

=   I    £[cos(wt)-l],    0<ojt<x,  (3.51) 

piston  1       * 

eCOS(ut),       TT<Olt. 

The  oscillation  frequency  is  f=l  Hz  (w=27rf=2T  rad/s).  The  geometric  parameters  are 
a=lxl0"3  m,  e=la,  and  L=3a.  The  nondimensional  parameters  are  7  =  1.4,  Pr=0.72, 
a=0.6559  and  M=1.846xl0~5.  It  is  known  from  thermodynamics  that  for  such  small 
dimensions  and  slow  oscillations  the  thermodynamic  quantities  should  be  uniform  within 
the  cylinder,  and  the  quantities  can  be  determined  analytically  by  the  adiabatic 
relationships: 
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Fig.  3.9   A  piston-cylinder  set-up 
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3.10   Computed  hydrodynamic  pressure  distributions  along 
the  x  axis  for  the  piston-cylinder  set-up 
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2l  =  (_£)T  =  (iLyr  .  (X)T-i  (  (3.52) 

Po" 


V  P„  T„ 


where  V  is  the  total  volume  and  the  subscript  "o"  denotes  initial  values. 

This  case  was  computed  numerically  on  a  machine  with  significant  digits  of  seven 
for  single  precision  and  of  16  for  double  precision.  It  was  first  computed  with  single 
precision  and  without  split  pressures  (denoted  as  Computation  No.  1),  next  it  was 
computed  with  double  precision  and  still  without  split  pressures  (Computation  No.  2), 
and  finally  the  same  case  was  computed  with  single  precision  and  with  split  pressures 
(Computation  No.  3).  The  numerical  pressure  distributions  along  the  x  axis  are  plotted 
against  each  other  in  Fig.  3.10.  The  pressure  shown  in  the  figure  is  actually  the 
pressure  difference  with  respect  to  the  point  (x,y)  =  (1.5a,0),  therefore  the  pressure  is  the 
hydrodynamic  pressure  ph.  As  indicated  by  the  figure,  at  time  o)t=ir/4,  the 
hydrodynamic  pressure  from  Computation  No.  1  oscillates  wildly  along  the  x  axis,  but 
at  o)t=2ir/3,  it  is  just  uniform.  None  of  these  two  distributions  is  physically  realistic. 
When  the  computation  accuracy  is  increased  by  using  double  precision  (Computation  No. 
2),  the  hydrodynamic  pressure  distribution  is  neither  oscillating  nor  uniform.  Such  an 
observation  reveals  that  for  problems  like  this,  the  pressure  field  cannot  be  correctly 
resolved  by  using  the  computational  algorithm  of  Computation  No.  1,  and  using  double 
precision  in  the  calculations  (Computation  No.  2)  is  indeed  one  remedy.  Figure  3.10 
clearly  shows  that  Computation  No.  3  yields  the  same  pressure  as  Computation  No.  2 
does.  Since  Computation  No.  3  uses  single  precision,  it  saves  computer  CPU  time  by 
about  30%,  memory  and  storage  by  almost  50%  over  Computation  No.  2.   This  proves 
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the  significant  benefits  given  by  using  the  pressure-splitting  algorithm,  and  verifies  the 

analysis  presented  in  the  last  chapter. 

The  instantaneous  thermodynamic  pressure,  p'(t),  in  the  third  oscillation  cycle  is 
shown  in  Fig.  3.11.  To  help  compare  the  magnitudes  of  ph  and  p1,  p'  in  this  figure  is 
normalized  with  the  same  quantity  as  that  for  ph.  By  comparing  Figs.  3.10  and  3.11,  it 
can  be  seen  that  the  magnitude  of  the  thermodynamic  pressure  is  greater  than  that  of  the 
hydrodynamic  pressure  by  nine  orders.  It  is  for  this  reason  that  Computation  No.  1  with 
only  seven  significant  digits  can  not  accurately  resolve  the  pressure  field  when  the 
thermodynamic  and  hydrodynamic  pressures  are  not  treated  separately  in  the  solution. 

In  Fig.  3.12,  the  numerical  temperature  distribution  along  the  x  axis  is  plotted 
against  the  analytical  solution.  At  a  given  time,  the  analytical  temperature  distribution 
is  uniform,  as  is  presented  by  the  horizontal  line  in  the  figure.  The  temperature  from 
Computation  No.  1  is  not  uniform,  but  varies  along  the  x  axis.  The  incorrect 
temperature  is  produced  by  the  erroneous  velocity  field  resulting  from  the  pressure  field 
that  is  not  accurately  resolved  numerically.  It  can  also  be  seen  that  the  numerical 
temperature  distribution  either  from  Computation  No.  2  or  from  Computation  No.  3  is 
in  perfect  agreement  with  the  analytical  solution,  which  again  shows  that  the  proposed 
pressure-splitting  technique  not  only  gives  the  right  solutions,  but  also  saves  computer 
CPU  time,  memory  and  storage  substantially. 

In  the  second  test  case,  the  set-up  is  similar  to  Model  Al  except  that  the  piston 
surfaces  are  isothermal.  Geometric  parameters  are  (refer  to  Fig.  1.3a):  a=lxl03  m, 
b=4a,  c=3a,  e=2a  and  L=6a.   The  wall  temperature  (including  the  piston's)  is  kept  at 


61 


3.5x10' 

\                                                                         /" 

i 

3.0x10" 
2.5x10s 

\                                                                /       - 

"a. 

2.0x10" 

\                                       / 

1.5x10' 

0  100  200  300 

cot  (degrees) 
Fig.  3.11    Thermodynamic  pressure  in  the  third  cycle  for  the  piston-cylinder  set-up 


x 

I- 


'      '      '      '      1 1      '      ' 

0,98 

<.)t=ir/4 

■ 

■ 

Symbol 

cot          Computation 

X 

n/4           No.  1 

0.96 

+ 

2ti/3         No.  1 

o 

nJ4           No.  2  or  3 

A 

2;i/3          No.  2  or  3 



Analytical 

U.94 

- 

No.1   Single  precision,  Without  pressure-splitting 

No.2  Double  precision,  Without  pressure-splitting 
No. 3  Single  precision,  Without  pressure-splitting 

0.92 

■ 

0.90 

Analytical;  No.2,  or  3 

art=2n/3                                       »♦♦••♦•♦« 

. 

,.♦*•♦♦•***,*  ****Vio. 

1 

1 
x/a 


Fig.  3.12   Temperature  distributions  along  the  x  axis  for  the  piston-cylinder  set-up 
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constant  Tw=288.15  K  (therefore  Tc=Th=Tw).  Pistons  oscillate  sinusoidally  with  the 
right  piston  leading  the  left  by  90°  in  phase  (yjo=90°)  with  an  oscillation  frequency  of 
f=2.32Hz.  Air  is  the  working  medium.  At  t=0,  p„=l  atm,  T0=  288.15  K,  p0=  1.225 
kg/m3.  The  corresponding  nondimensional  parameters  are  7  =  1.4,  Pr=0.72,  a=l  and 
M=4.29xl0'5. 

In  this  calculation,  the  pressure-splitting  algorithm  is  used.  For  the  computational 
domain,  which  is  only  the  upper  half  geometry,  a  91x31  grid  system  is  used.  The  grid 
distribution  is:  31x31  in  the  left  chamber,  31x11  in  the  channel,  and  31x31  in  the  right 
chamber.  The  time  step  is  chosen  such  that  each  cycle  is  discretized  into  360  uniform 
intervals,  i.e.,  At=  l/(360f  )=  1.2xl03  sec.  It  takes  about  four  cycles  for  the  solution  to 
reach  the  time-periodic  state. 

The  P-V  diagram  for  this  calculation  is  shown  in  Fig.  3.13.  This  figure  shows 
that  the  minimum  pressure  is  pn,i„/p„= 0.7099  and  the  maximum  is  p,^, /p„=  1.6420, 
which  can  be  checked  by  examining  two  extreme  conditions,  namely  isothermal  and 
adiabatic.    For  an  isothermal  condition, 

p/p„  =  V0/V  , 
so,  Pmi„/p„  =  V0/Vm„  =  0.7262,  and  p^/p,,  =  v^V^  =  1.6055. 

For  an  adiabatic  condition, 

p/p„  =  (V„/V)\ 
so,  p^/p,,  =  (v^v^1  =  0.6389,  and  p„„x/p0  =  (V^V^)7  =  1.9402. 

The  numerical  maximum  should  be  larger  than  the  adiabatic  counterpart  and  less  than  the 
isothermal  one,  whereas  the  numerical  minimum  should  be  larger  than  the  isothermal  one 
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Fig.  3.13   P-V  diagram  of  the  second  test  case 
and  less  than  the  adiabatic  one;  moreover,  numerical  values  should  be  closer  to  the 
isothermal  values  than  to  the  adiabatic  ones  because  of  the  boundary  condition  of 
isothermal  walls.    That  is  indeed  the  case: 

minimum:        0.6389  (adiabatic)  <  0.7099  (numerical)  <  0.7262  (isothermal), 

maximum:       1.6035  (isothermal)  <  1.6420  (numerical)  <  1.9402  (adiabatic). 

The  correctness  of  the  numerical  solution  is  also  supported  by  checking  the 

overall  energy  balance.    First  define  the  power  input— the  mean  rate  at  which  work  is 

done  on  the  gas  by  some  external  force— as 
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I      ft   f-„HVI* 


P  =  -   |{/pdV)dt, 


0      V(t) 

where  t  is  the  period  of  piston  oscillations,  i.e.  t  =  2ir/oj. 

The  mean  rate  at  which  heat  is  transferred  out  of  some  domain  is 


T 

q  =  I  j  {  j  -kVT-ndS)dt  ,  <3-54) 


where  Swe,(t)  represents  the  solid  surface  that  is  in  contact  with  the  gas,  and  the  surface 
may  change  with  time. 

The  rate  at  which  the  gas  "absorbs"  heat  through  the  left  chamber  walls  plus  the 
piston  surface  is  Q,efl=2.4589D  watts  (here  D  is  the  depth  in  the  third  dimension  in 
meters);  the  rate  at  which  heat  is  rejected  out  of  the  gas  through  the  channel  walls  is 
Qohlnoel=0.5516D  watts,  and  that  through  the  right  chamber  walls  plus  the  piston  surface 
is  Qright=2.7297D  watts.  The  power  input  is  P=0.8415D  watts.  According  to  the  first 
law  of  thermodynamics,  the  following  expression  holds, 

Qlefl  "  Vchannel  "   Vrighl   *"   ""■ 

Then,   LHS  =  2.4589  -  0.5516  -  2.7297  =  -0.8224, 

RHS  =  -0.8415. 
Therefore,  there  is  only  a  very  small  relative  error  of  2.3%  in  the  total  energy  balance. 

The  preceding  two  test  cases  have  verified  that  the  pressure-splitting  algorithm  is 
computationally  accurate  and  efficient,  and  in  the  mean  time  they  have  also  validated  the 
computer  code. 


CHAPTER  4 
RESULTS  AND  DISCUSSION 


By  using  this  code,  investigations  were  made  on  the  models'  cooling  capacity  and 
efficiency  in  terms  of  working  gases,  oscillation  frequencies,  channel  length,  ratios  of 
the  chamber  width  to  the  channel  width,  etc.  The  velocity  and  temperature  fields  were 
also  examined.  Discussion  is  made  on  the  criterion  for  the  model  to  result  in  good 
regeneration.  Table  4. 1  lists  the  cases  that  were  studied.  Note  that  to  make  the  number 
of  parameters  appearing  in  the  table  as  small  as  possible,  yet  adequate  to  present  each 
case,  those  parameters  that  can  be  obtained  directly,  or  via  auxiliary  relations  if  needed, 
are  omitted  there.  For  instance,  parameters  that  do  not  need  such  relations  are:  the  ratio 
of  specific  heats  (7=1.4  for  air  and  7  =  1.7  for  helium),  Prandtl  number  (Pr=0.72  for 
air  and  Pr=0.67  for  helium),  and  the  constant  T,  used  in  the  Southerland's  formula  (Eq. 
(3.9))  (T„=  110.4  K  for  air  and  T„=61.67  K  for  helium).  One  of  the  parameters  that 
need  such  relations  is  the  reference  density  p, ,  which  can  be  obtained  via  the  equation 
of  state  as  Pr=pr/(RTr)  with  R  being  287  J/(Kg-K)  for  air  and  2077.2  J/(Kg-K)  for 
helium.  For  the  meanings  of  the  geometric  parameters  listed  in  Table  4.1,  please  refer 
to  Figs.  1.3  and  1.4. 
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4.1  Velocity  and  Thermodynamic  Fields  of  a  Typical  Case 
In  this  section,  the  computations  and  results  of  Case  I  (refer  to  Table  4.1)  are 
examined.  Computations  were  made  on  the  Cray  C90  supercomputer  of  the  Pittsburgh 
Supercomputer  Center.  Different  grids  and  time  step  sizes  were  first  tested  to  make  sure 
that  they  were  fine  enough  for  the  numerical  solutions  to  be  accurate  enough.  It  was 
identified  that  for  this  particular  case  a  reasonable  choice  was:  a  grid  system  of  2153 
points  (31x31  in  the  left  chamber,  21x1 1  in  the  channel,  and  31x31  in  the  right  channel) 
and  a  time  step  size  of  At  =  r/360  ( =  2.78x10 3  sec),  for  further  refining  either  the  grid 
or  the  time  step  size  did  not  result  in  any  noticeable  changes  in  the  numerical  solutions. 
The  grids  at  o)t  =  ir/4  and  ut=ir  for  this  case  can  also  be  found  in  Figs.  3.1  (a)  and  (b) 
respectively,  though  these  figures  originally  served  another  purpose  there.  Solution 
iteration  numbers  used  were:  8  iterations  for  solving  the  u-momentum  equation,  8  for  the 
v-momentum  equation,  40  for  the  p'  equation  and  8  for  the  energy  equation;  and  300 
repetitions  of  the  above  procedures  at  a  give  time  level.  Under  these  given  conditions, 
it  took  about  5472  second  CPU  time  to  run  one  cycle  on  the  Cray  C90  machine. 

Figures  4. 1  (a)  and  (b)  plot  the  transient  processes  of  the  thermodynamic  pressure 
and  the  pressure  drop  between  the  left  and  right  pistons.  These  figures  reflect  the  slow- 
start  velocity  and  temperature  conditions  in  the  initial  quarter  cycle  discussed  at  the  end 
of  the  last  chapter.  Figure  4. 1  (b)  shows  that  the  (hydrodynamic)  pressure  drop  adjusts 
itself  to  time-periodic  state  only  in  few  cycles.  The  thermodynamic  pressure  oscillates 
in  the  range  of  0.5  to  2.5  atms,  whereas  the  hydrodynamic  pressure  drop  is  only  in  the 
range  of  -2xl0"6  to  2xl0"6  atms,  with  the  latter  being  only  about  one  millionth  of  the 
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(a)  Thermodynamic  pressure 

a=1  mm,  b=8a,  c=8a,  e=7a,  L=1 2a,  f=1  Hz,  cp0=90° 
Helium,  p=1  atm,  T  =288.1 5K,  Tt=0.9T, ,  T>1 .1T, ,  a=0.229 

2.5 


t  /  x   cycle 


(b)  Pressure  drop 

a=1mm,  b=8a,  c=8a,  e=7a,  L=12a,  f=1Hz,  <po=90° 
E  Helium,  p=1atm,T>288.15K,T>0.9Tr,T,,=1.1T,,a=0.229 

*         4x1 0E 


Q. 


-2x10 


Fig.  4.1    Transient  processes  of  pressures, 
(a)  Thermodynamic  pressure;    (b)  Pressure  drop  between  pistons. 
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former. 

A  way  to  check  if  numerical  solutions  have  reached  the  time-periodic  state  is  to 
compare  the  current  solutions  with  those  at  the  time  one  cycle  earlier.  The  quantity  a 
defined  below  is  such  a  measure. 

o,  -  -J .  (4.1, 

E^ifljl  +  i*STI» 

where  4>  is  a  general  notation  representing  any  specific  dependent  variable  of  u,  v,  p\ 
T  or  p.  Such  a  measure  is  objective,  because  it  takes  the  overall  domain  into  account 
by  the  summation  and  meanwhile  is  a  relative  measure  by  adding  the  denominator  in  the 
definition.  Figure  4.2  plots  the  quantities  a„,  ay ,  erph  and  aT  in  the  logarithmic  scale 
vs.  time.  At  any  given  time,  they  are  of  the  same  order  and  they  all  decrease  linearly 
as  time  increases.   The  lines  of  logo's  have  exactly  the  same  slope. 

Figures  4.3  (a),  (b),  (c)  and  (d)  give  the  velocity  vectors  at  time  instants  o>t=0, 
ir/2,  x  and  3x/2,  respectively.  Note  that  to  have  a  better  view  of  the  field  quantity  in 
the  configuration,  the  domain  shown  in  the  figures  is  actually  as  large  as  six  times  the 
computational  domain.  No  significant  circulation  region  is  observed  for  this  case  because 
of  the  small  Womersley  number  used. 

The  corresponding  hydrodynamic  pressure  contours  are  plotted  in  Figs.  4.4  (a), 
(b),  (c)  and  (d).  These  figures  show  that  when  neither  piston  is  extremely  close  to  the 
channel  exit,  pressure  drops  mainly  occur  in  the  channels,  otherwise  the  pressure  drop 
in  the  transverse  direction  within  the  chamber  is  comparable  to  that  in  the  channel.  The 
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a=1  mm,  b=8a,  c=8a,  e=7a,  L=1 2a,  f=1  Hz,  cp0=90° 
Helium,  p  =1  atm,  T=288.1 5K,  T0=0.9Tr ,  T>1 .1T, ,  a=0.229 


Fig.   4.2   Convergence  history  to  time-periodic  state 
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total  pressure  drop  in  space  is  only  a  few  millionths  of  an  atmosphere  at  most.    It  will 
be  seen  later  that  this  is  in  contrast  with  the  thermodynamic  pressure,  which  changes  with 
time  by  an  amount  in  the  order  of  one  atmosphere. 

Figures  4.5  (a),  (b),  (c)  and  (d)  demonstrate  the  temperature  contours.  At  all  the 
times,  the  temperature  of  the  gas  within  the  channel  is  very  close  to  the  local  channel 
wall  temperature,  indicating  that  the  channel  functions  as  a  regenerator  properly.  At 
ut=0  (Fig.  4.5  (a)),  which  is  at  the  compression  stage,  the  temperature  in  the  right 
chamber  is  much  higher  than  that  of  the  right  chamber  wall  (1.1  nondimensional),  and 
a  value  as  high  as  1.625  is  observed.  Inversely,  at  wt  =  x  (Fig.  4.5  (c)),  which  is  at  the 
expansion  stage,  the  temperature  in  the  left  chamber  is  much  lower  than  that  of  the  left 
chamber  wall  (0.9),  and  a  value  of  0.625  is  observed  there.  On  average  over  time  and 
space  the  mean  gas  temperature  in  the  right  chamber  is  higher  than  that  of  the  chamber 
wall,  so  heat  is  rejected  out  of  the  gas  through  the  right  chamber  wall;  on  the  other 
hand,  the  mean  gas  temperature  in  the  left  chamber  is  lower  than  that  of  the  chamber 
wall,  and  heat  is  absorbed  by  the  gas  there. 

Density  contours,  shown  in  Fig.  4.6  (a),  (b),  (c)  and  (d),  indicate  that  density 
changes  significantly  with  time  and  space.  For  the  time  instants  chosen,  the 
nondimensional  density  takes  values  from  0.55  to  2.4.  At  any  time,  the  patterns  of 
density  contours  look  like  those  of  temperature  contours  due  to  the  reciprocal  relationship 
between  density  and  temperature  at  a  given  time  (Eq.  (3.9)). 

Figure  4.7  plots  the  P-V  diagram  for  this  case.  The  direction  of  the  loop  is 
counterclockwise,  indicating  that  some  external  force  does  work  on  the  gas.   The  power 
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Fig.  4.7  P-V  diagram  of  Case  1 
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input  is  P  =  14.075D  watts  (as  mentioned  before,  the  symbol  "D"  represents  the  third 

dimension  in  meters).  The  rate  at  which  the  gas  absorbs  heat  through  the  solid  surfaces 
having  cold  temperature  (Tc)  is  defined  as  the  cooling  capacity  and  is  denoted  by  Qc, 
while  the  rate  at  which  heat  is  rejected  from  the  gas  through  the  solid  surfaces  having  hot 
temperature  (Th )  is  defined  as  the  heating  rate  and  denoted  by  Qh.  Heat  can  also  be 
conducted  into  or  out  of  the  channel  walls  along  which  the  temperature  distribution  is 
linear.  This  is  an  undesirable  leakage,  the  rate  of  which  is  denoted  by  Qletk.  The 
quantities  here  are:  Qc  =  11.574D  watts,  Qh  =  21.643D  watts,  and  Qlalt  =  3.990D 
watts.    So  the  total  rate  of  energy  into  the  system  is 

p  +  Qc  =  14.075D  +  11.574D  =  25.649D, 
and  the  total  rate  of  energy  out  of  the  system  is 

Qh  +  Qfc„t  =  21.643D  +  3.990D  =  25.633D. 
Hence  the  overall  energy  balance  is  excellent. 

In  the  above  discussions,  the  factor  "D"  and  unit  "watt"  always  go  with  the 
quantities:  CL ,  Qh  ,Qleak  and  P.  The  factor  and  units  will  be  omitted  in  the  following  for 
simplicity. 

4.2  Model  Comparison 

To  compare  Models  Al  and  A2,  consider  Cases  2  and  3.   As  listed  in  Table  4. 1, 

all  the  conditions  for  these  cases  are  the  same  except  that  in  Case  2  the  temperature  along 

the  entire  channel  wall  is  linear,  while  in  Case  3  the  channel  wall  temperature  is  linear 

only  along  the  middle  segment  of  the  channel  and  the  wall  temperature  is  constant  along 
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either  end  segment.  Table  4.2 
compares  the  performances  of  the  two 
cases.  The  negative  value  of  Q,^  in 
Case  3  means  that  heat  is  added  into 
the  gas  through  the  channel  walls 
having  a  linear  temperature 
distribution.  In  the  table,  COP  stands 
for  the  coefficient  of  performance, 
which  is  defined  as  (Walker  1973) 


Table  4.2  Comparison  of  Models 
Al  and  A2 


Case  2 
Model  Al 

Case  3 
Model  A2 

Q. 

12.7403 

11.8207 

Qh 

17.4084 

19.3578 

Vleak 

2.4482 

-0.3929 

P 

7.1340 

7.1626 

COP 

1.7859 

1.6503 

V, 

39.69% 

36.67% 

(4.2) 


COP  =  Cooling  capacity  /  Power  input  =  Qc/  P. 
The  COP  of  Case  2  is  1.7859,  and  that  of  Case  3  is  1.6503;  the  latter  is  only  slightly 
lower  than  the  former.  According  to  the  Carnot  theorem,  of  all  the  thermodynamic 
cycles  with  common  values  for  the  maximum  and  minimum  temperatures,  pressures  and 
volumes,  the  Carnot  cycle  has  the  maximum  COP.  For  a  cooling  device,  the  COP  of 
a  Carnot  cycle  is  (Walker  1973) 

(COP)™  =  Tc  /  (T,  -  Tc ).  (4.3) 

In  Table  4.2,  the  symbol  i)r  represents  the  relative  efficiency  defined  by  (Walker  1973) 

n,  =  Actual  COP  /  Carnot  COP  =  COP  /  (COP)™  .  (4.4) 

The  relative  efficiency  of  an  actual  device  usually  varies  between  40%  and  a  few  percent 
depending  on  the  size  of  the  device  and  the  operation  conditions.  Here  the  relative 
efficiency  of  either  Case  3  or  4  is  only  slightly  below  40%.  In  terms  of  the  cooling 
capacity,  power  input  and  hence  COP,  Models  Al  and  A2  have  almost  the  same 
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performance  under  the  given  conditions. 

Next  compare  Models  A2  and  B2.   Case  3  uses  the  single-channel  model  A2  and 

Case  4  uses  the  multichannel  model  B2,  with  all  other  conditions  being  the  same.   It  can 

be  seen  from  Table  4.3  that  Case  4  gives  a  cooling  capacity  of  11.3202,  which  is  quite 

close   to   Case   3's    11.8207.      Case   4 

Table  4.3  Comparison  of  Models 
requires  a  power  input  of  1 1.3809,  while  A2  and  B2 

Case  3  requires  only  7. 1626.   As  a  result. 

Case  4  has  a  COP  as  low  as  0.9947.  The 

use  of  the  multichannel  model  in  Case  4 

makes   the   compression  and   expansion 

processes  closer  to  adiabatic  than  in  Case 

3  where  single-channel  model  is  used,  so 

Case  4  should  have  a  larger  range  of 

pressure    variations    despite    the    same 

volumetric  changes  for  the  two  cases.     Computed  data  (not  tabulated  here)  indicate 

pressure  limits  of  [0.64,  2.12]  (atm)  for  Case  3  and  limits  of  [0.57,  2.22]  for  Case  4. 

That  explains  why  the  latter  case  needs  a  larger  power  input  than  the  former.    The 

compression  and  expansion  processes  in  Case  4  produce  lower  temperature  in  the  left 

chamber  and  higher  temperature  in  the  right  chamber,  respectively,  than  in  Case  3;  the 

cooling  or  heating  of  Case  4,  however,  is  still  close  to  the  counterpart  of  Case  3,  because 

the  heat  exchanging  area  of  Case  4  is  much  smaller  than  that  of  Case  3.    Even  though 


Case  3 
Model  A2 

Case  4 
Model  B2 

Qc 

11.8207 

11.3202 

Q„ 

19.3578 

22.9956 

Qleak 

-0.3929 

-0.3421 

P 

7. 1626 

11.3809 

COP 

1.6503 

0.9947 

V, 

36.67% 

22.10% 

85 
the  multichannel  model  has  a  lower  COP  than  the  single-channel  model  does,  the  former 

is  more  practical  from  an  engineering  viewpoint. 


4.3  Comparison  of  Using  Helium  and  Air 

The  most  widely  used  working  gases  for  Stirling-type  machines  are  air,  helium 

and  hydrogen.   The  selection  of  a  gas  for 

Table  4.4  Comparison  of  helium  and  air 
a  specific  Stirling  cooling  device  not  only 

depends      on      the      thermodynamic 

performance  of  the  gas,  but  also  depends 

on  other  factors  such  as  cost,  availability, 

safety  and  so  on  (Reader  and   Hooper 

1983).     With  the  computer  code,  it  is 

very      convenient      to      compare      the 

performances  of  working  gases.  Consider 

Cases  4  and  5,  in  which  geometric  parameters,  temperature  boundary  conditions,  piston 

oscillation  frequencies  and  amplitudes  are  the  same,  but  helium  is  the  working  gas  in  the 

former  and  air  in  the  latter.    From  Table  4.4  it  can  be  seen  that  helium  produces  a 

cooling  capacity  twice  as  large  as  air  does,  while  the  power  input  needed  for  one  case 

is  close  to  the  other.    As  a  result,  helium  is  about  twice  as  efficient  as  air  under  the 

given  conditions. 

There  are  mainly  two  reasons  contributing  to  this.     First,  the  monatomic  gas 

helium  has  a  ratio  of  specific  heats  of  1.67,  which  is  higher  than  air's  1.40.   So,  for  the 


Case  4 
Helium 

Case  5 
Air 

Qc 

11.3202 

5.2169 

Qh 

22.9956 

14.1773 

Qick 

-0.3421 

1.2477 

p 

11.3809 

10.2219 

COP 

0.9947 

0.5104 

V, 

22.10% 

11.34% 
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same  amount  of  expansion,  more  cooling  is  produced  by  helium  than  by  air.  Second, 
helium  has  a  smaller  Prandtl  number  than  air  (0.67  versus  0.72)  and  a  kinematic 
viscosity  eight  times  as  large  as  air.  As  result,  the  channel's  heat  regeneration  for 
helium  is  better  than  for  air.  And  a  detailed  discussion  on  the  regeneration  will  be  found 
in  Section  4.8. 


4.4  Phase  Angle  Effect 
In  order  to  examine  the  effect  of  phase  angle  on  the  performances  of  the  models, 
five  runs  (Cases  6  through  10  in  Table  4.1)  were  made  with  phase  angles  of  <?„  =  45°, 
75°,  90°,  105°  and  135°,  respectively,  with  other  parameters  being  maintained  constant. 
Figure  4.8  plots  the  cooling  capacity  (QJ,  rate  of  heat  rejected  (Q,),  power  input  (P) 
and  COP  versus  the  phase  angle  <p0 .  The  figure  shows  that  each  of  the  four  quantities 
takes  its  maximum  in  the  range  of  phase  angles  covered,  and  that  the  phase  angle 
corresponding  to  the  maximum  value  of  one  quantity  may  differ  from  that  of  another. 
For  instance,  the  power  input  has  a  maximum  at  <p0  —  110°,  the  cooling  capacity  a 
maximum  at  <p„  -  96°,  and  COP  a  maximum  at  <p0  —  75°.  The  optimum  phase  angle 
for  the  cooling  power  is  not  optimum  for  COP.  Generally  speaking,  the  phase  angles 
corresponding  to  the  maxima  are  not  too  far  off  90°,  and  each  of  the  four  quantities 
shown  in  the  figure  is  not  very  sensitive  to  changes  in  ip0  near  the  optimum  phase  angle. 
These  observations  based  on  the  numerical  simulations  are  in  agreement  with 
experimental  data  (Reader  and  Hooper  1983,  pp.  75-77). 
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a=1  mm,  b=6a,  c=3a,  e=2a,  L=1  Oa,  f=1  Hz 
Helium,  p=1atm,  T  =288.1 5K,  T.-09T,,  T>1 .1T,,  a=0.229 


75  90  1 05 

Phase  angle  q>0  (°) 
Fig.  4.8  Phase  angle  effects 


4.5  Impact  of  Varying  Refrigeration  Temperatures  on  Performances 
To  understand  the  effects  of  changing  the  refrigeration  temperature  on  the 
performances,  several  runs  (Cases  2  and  11  through  15)  were  made,  in  which  the 
temperature  difference  between  the  cold  and  hot  chambers  was  changed  to  a  large  extent 
while  all  other  parameters  were  kept  unchanged.  In  Fig.  4.9,  the  cooling  capacity, 
power  input  and  COP  are  plotted  versus  the  ratio  of  the  cold  to  hot-end  temperature, 
Tc/Th .  Also  plotted  in  the  figure  is  the  mean  rate  at  which  heat  is  conducted  and 
convected  into  the  cold  chamber  through  the  gas  at  the  cross  section  of  the  channel  exit. 
This  rate  is  regarded  as  a  loss  (denoted  by  Qi„s),  because  it  will  consume  some  cooling 
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a=1  mm,  b=6a,  c=3a,  e=2a,  L=1  Oa,  f=1  Hz,  if>„=90° 
Helium,  p=1atm,a=0.229 


Fig.  4.9  Refrigeration  temperature  effects 


made  by  expansion. 

The  cooling  capacity,  power  input,  the  loss  and  the  COP  all  vary  monotonously 
with  the  temperature  ratio.  When  the  refrigeration  temperature  rises,  the  cooling 
capacity  increases,  but  the  power  input  and  the  loss  decrease,  therefore  the  COP 
increases,  too.  The  closer  the  cold  chamber  temperature  is  to  the  hot-end  temperature 
(usually  the  ambient),  the  better  is  the  performance.  For  instance,  when  Tc/Th  equals 
0.76,  the  COP  is  about  1.  However,  as  Tc/Th  is  increased  to  0.9,  the  COP  is  close  to 
4.  Such  phenomena  were  also  reported  by  McDonald  and  coworkers  (McDonald  et  al. 
1994)  based  on  their  design  analysis  and  experimental  data  with  an  actual  Stirling 
refrigerator. 
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When  the  cold  chamber  wall  temperature  is  set  too  low,  cooling  is  no  longer 

possible.    Instead,  heat  is  rejected  out  of  the  gas  in  the  left  chamber  through  the  walls. 

That  is  indicted  by  the  negative  cooling  capacity  in  the  figure.  If  the  hot-end  temperature 

Th  is  300  K,  the  lowest  obtainable  refrigeration  temperature  will  be  about  189  K  under 

the  operation  conditions  considered.  The  large  heat  loss  at  low  refrigeration  temperatures 

can  consume  all  the  cooling  produced  by  expansions. 


4.6  Changing  Channel  Length 
Increasing  the  channel  length  (or  the  regenerator  length)  will  decrease  the 
temperature  gradient  along  the  channel,  and  therefore  can  reduce  the  conduction  loss. 
Comparing  Case  2  and  Case  16  (Table 


Table  4.5  Comparison  of  cases  with 
different  channel  lengths 


4.5)  will  demonstrate  that.    The  channel 

length  in  the  former  is  24a  and  that  in 

the  latter  is  12a,  and  other  parameters  in 

both  cases  are  the  same.    The  loss  with 

the  longer  channel  is  0.6131,  which  is 

less  than  1 .2443-the  counterpart  with  the 

shorter  channel.    However,  the  cooling 

capacity    with    the    longer    channel    is 

12.7403,  which  is  surprisingly  even  less  than  the  value  (13.745)  of  the  cooling  capacity 

with  the  shorter  channel.    The  reason  for  this  is  that  although  the  loss  is  reduced  by 

lengthening  the  channel,  meanwhile  the  volume  that  cannot  be  swept  by  the  pistons-the 


Case  2 
L=24a 

Case  16 
L=12a 

Vloss 

0.6131 

1.2443 

Qc 

12.7403 

13.7450 

P 

7.1340 

8.6155 

COP 

1.7859 

1.6000 
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dead  volume-is  also  increased,  therefore  the  cooling  produced  by  expansion  decreases. 
The  net  effect  (i.e.,  the  available  cooling  Qc)  depends  on  these  two  competing  factors 
under  specific  operation  conditions.  Despite  such  complexity,  it  is  always  true  that 
increasing  the  channel  length  also  increases  the  COP  because  the  loss  is  reduced.  Table 
4.5  does  show  that  the  COP  with  the  longer  channel  is  1.7859  but  that  with  the  shorter 
one  is  only  1.6. 

4.7  Effects  of  Changing  the  Ratio  of  Chamber  Width  to  Channel  Width 
Under  given  conditions,  if  the  ratio  of  chamber  width  (2b)  to  channel  width  (2a), 
i.e.  b/a,  is  too  small,  the  gas  in  the  cold  chamber  is  not  effectively  separated  from  that 
in  the  hot;  besides,  the  effective  cooling  surface  is  also  too  small.  On  the  other  hand, 
if  the  ratio  is  too  large,  the  gas  in  one  chamber  cannot  be  promptly  transferred  to  the 
other  when  needed.  Either  situation  may  degrade  the  cooling  performance.  For  the 
multichannel  models,  in  addition  to  the  COP,  another  appropriate  measure  of  the 
performance  is  the  cooling  capacity  per  unit  area  of  the  chamber  cross  section,  i.e.  Qc 
/(2bD),  which  indicates  how  economically  the  device  uses  space.  Four  cases  were  run 
with  the  following  values  of  b/a:  4  (Case  18),  8  (Case  1),  12  (Case  18)  and  15  (Case 
19).  Figure  4. 10  plots  COP  and  Qc  /(2bD)  versus  b/a.  The  figure  shows  that  both  COP 
and  Qc/(2bD)  vary  with  b/a,  and  there  exists  an  optimum  value  of  b/a  for  either  quantity. 
When  b/a  is  approximately  6.4,  the  efficiency  of  performance  takes  its  maximum  value 
of  about  0.84.  When  b/a  is  about  1 1.5,  the  quantity  Qc/(2bD)  reaches  a  maximum  value 
of  about  750  Watts/m2 . 
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a=1  mm,  c=8a,  e=7a,  L=12a,  <po=90°,  f=1  Hz 
Helium,  p=1atm,  T=288.15K,  T>0.9Tr,  T„-1 .1T„  <x=0.229 
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Fig.  4. 10  Effects  of  changing  the  ratio  of  chamber  width  to  channel  width 

4.8  Discussion  on  Regeneration 
It  is  known  that  the  regeneration  efficiency  can  greatly  affect  the  performances 
of  a  Stirling  device.  This  section  discusses  the  condition  required  for  the  channel  with 
a  linear  temperature  distribution  to  provide  adequate  regeneration.  For  good 
regeneration,  there  should  be  enough  time  for  the  gas  flowing  through  the  channel  to 
adjust  its  temperature  to  the  that  of  the  local  channel  wall.  The  time  scale  of  the  thermal 
diffusion  across  the  channel  is 
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t     _a2  _  a^Pr  (45) 

ldif  ' 

K  V 

where  k  is  the  thermal  diffusivity.  The  residence  time  during  which  the  gas  passes 
through  the  channel  and  exchanges  heat  with  the  channel  is 

t     ■  —  ,  (4.6) 

res  v       > 

in  which  U„,  is  the  mean  magnitude  of  the  gas  velocity  within  the  channel.  This  quantity 
is  estimated  by  using  the  ratio  of  the  chamber  width  (2b)  to  the  channel  width  (2a)  and 
the  mean  magnitude  of  the  piston  velocity  Upm.  Since  the  piston  velocity  is  in  the  form 
of  fi«cos(cdt+<p),  the  quantity  Upm  is  ea/y[2.   Then  one  has 

Um  -  ^X  -  *  «•  (4.7) 

a    m       a  ft 

It  is  necessary  for  tres  to  be  greater  than  tdit,  i.e., 

Substituting  Eqs.  (4.5)-(4.7)  into  Eq.  (4.8)  and  then  rearranging  terms  yield 

P  =  J_^la2Pr  <  1  .  (4.9) 

v^aL 

The  dimensionless  parameter  0  defined  in  Eq.  (4.9)  is  a  combination  of  two  other 

parameters  related  to  the  geometry  (i.e.  b/a  and  e/L),  flow  oscillation  (i.e.  a)  and  heat 

transfer  (i.e.  Pr).    For  effective  regeneration,  0  should  be  kept  far  below  unity. 

Compare  Cases  8  and  20.   The  parameter  0  in  the  former  is  0.0298,  and  that  in 

the  latter  is  0.298.   Such  an  increase  in  /3  is  made  possible  by  increasing  the  oscillation 
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frequency  by  nine  times.  Figures  4.11  (a)  and  (b)  plot  the  temperature  distributions 
along  the  x  axis  at  eight  time  instants  for  the  two  cases  respectively.  For  the  case  with 
the  smaller  0,  at  all  the  time  instants  and  within  almost  the  entire  channel,  the  gas 
temperature  is  always  very  close  to  the  local  channel  wall  temperature.  This  indicates 
that  the  channel  functions  well  as  a  regenerator.  However,  for  the  case  with  the  larger 
0,  the  temperature  of  the  gas  within  the  channel  deviates  significantly  from  that  of  the 
local  wall.  As  a  result,  the  case  with  (3=0.0298  has  a  COP  of  1.8606,  but  the  case  with 
(3=0.298  has  merely  a  COP  of  0.9265.  The  comparison  made  here  clearly  shows  the 
role  played  by  the  parameter  (3  and  the  need  to  keep  this  well  below  unity  in  a  typical 
channel  regenerator. 
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(a) 

a=1  mm,  b=6a,  c=3a,  E=2a,  L=1  Oa,  (po=90°,  f-1  Hz,  a=0.229 
Helium,  p=1atm,  T=288.15K,  Tc=0.9Tr, Th=1 .1T, 
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a=1  mm,  b=6a,  c=3a,  E=2a,  L=1  Oa,  <p0=90°,  f=1 0Hz,  a=0.724 
(b)  Helium,  p=1atm,Tr =288.1 5K,Tc=0.9Tr,Th=1.1Tf 


Fig.  4.11  Temperature  distributions  along  x  axis  for  two  different  /3's 
(a)  0  =  0.0298;    (b)  /3  =  0.298. 


CHAPTER  5 
CONCLUDING  REMARKS 


A  numerical  technique  was  devised  to  separate  the  hydrodynamic  and 
thermodynamic  pressures  for  low-Mach  number  variable-density  flows  with  heat  transfer, 
which  occur  in  reciprocating  thermodynamic  machines.  Separation  of  these  two 
pressures  different  in  orders  of  magnitudes  avoids  the  catastrophic  cancellations  in 
numerically  evaluating  the  pressure  gradients,  and  so  the  technique  can  increase 
computational  accuracy  and  efficiency  substantially  (potential  saving  on  CPU  time  is 
about  30%  and  that  on  storage  and  memory  is  50%).  It  can  be  easily  incorporated  to  any 
existing  pressure-based  numerical  algorithms. 

Two-dimensional  models  of  three-component  Stirling  microrefrigerators  were 
developed  in  this  study.  The  models  consist  of  all  three  heat  exchangers:  heater,  cooler 
and  regenerator.  The  heater  and  cooler  are  ended  with  pistons  oscillating  out  of  phase, 
and  connected  by  narrow  channels  with  a  linear  wall-temperature  distribution,  simulating 
the  regenerator.  Numerical  computations  were  made  for  the  time-periodic,  variable- 
density  flow  and  heat  transfer  within  domains  with  moving  boundaries  by  using  the 
SIMPLE  algorithm  incorporated  with  a  moving  grid  technique  and  with  the  newly 
devised  pressure-splitting  method.  Such  a  solution  methodology  was  proven  to  be 
feasible-accurate  and  efficient-for  the  complex  flow  and  heat  transfer  problems  in  the 
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Stirling-type  devices.  This  is  the  first  reported  attempt  to  systematically  discuss  the 
numerical  aspects  pertaining  to  the  simulations  of  Stirling  devices.  Results  based  on  the 
microrefrigerator  models  agreed  very  well  with  the  experimental  data  obtained  from 
actual  devices. 

The  model's  regeneration  and  cooling  performance  are  sensitive  to  the 
dimensionless  parameter  that  is  the  product  of  Womersley  number  squared,  Prandtl 
number,  the  relative  piston  strokes,  and  the  ratio  of  chamber  width  to  channel  width. 
This  parameter  should  be  kept  far  below  unity  for  the  connecting  channels  to  function 
well  as  a  regenerator. 

As  the  cold  chamber  wall  temperature  decreases,  the  cooling  capacity  decreases 
linearly,  but  the  power  input  increases  linearly,  and  consequently  the  performance 
degrades.  For  a  set  of  given  operation  conditions,  there  is  a  lowest  temperature  below 
which  no  cooling  is  possible.  This  is  caused  by  the  increasing  cooling  loss  through  the 
channel  as  the  cold  chamber  temperature  decreases.  Lengthening  the  channels  decreases 
such  loss,  but  meanwhile  increases  the  dead  space,  which  in  turn  reduces  the  cooling 
made  by  the  expansion.  The  net  effect  depends  on  these  two  competing  factors  under 
specific  conditions. 

The  phase  angle  affects  the  cooling  capacity  and  efficiency.  The  phase  angle  that 
results  in  the  maximum  cooling  capacity  may  not  necessarily  result  in  the  maximum 
efficiency.  The  optimum  phase  angle  for  either  cooling  capacity  or  efficiency  is  near 
90°. 

The  ratio  of  the  chamber  width  to  the  channel  width  also  affects  the  cooling 
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performance.     Under  given  operation  conditions,  there  is  an  optimum  ratio  for  the 

cooling  capacity  and  efficiency  respectively. 

Future  work  on  this  topic  should  be  in  the  following  directions.  First,  turbulence 
models  for  oscillatory  internal  flows  should  be  incorporated  into  the  computations  so  that 
simulations  can  also  be  made  for  larger  devices  and  for  higher  frequencies.  Second, 
model  the  regenerator  as  porous  medium  to  simulate  the  conjugate  heat  transfer  more 
accurately.  Third,  use  a  five-component  model  instead  of  the  three-component  model. 
Finally,  the  dimensions  of  the  model  should  be  increased  from  two  to  three.  These 
extensions  will  need  a  very  large  increase  in  computing  resources  of  both  computer 
memory  and  CPU  time,  and  it  will  be  necessary  to  use  parallel  computations  on  Cray 
supercomputers. 
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